##// END OF EJS Templates
schain
RSS
First Draft HDF5 IO module
Julio Valdez -
r514:f095b959308c
parent child
Show More
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (656 lines changed) Show them Hide them
@@ -0,0 +1,656
1 import numpy
2 import time
3 import os
4 import h5py
5 import re
6
7 from model.data.jrodata import *
8 from model.proc.jroproc_base import ProcessingUnit, Operation
9 from model.io.jroIO_base import *
10
11
12 class HDF5Reader(ProcessingUnit):
13
14 ext = ".hdf5"
15
16 optchar = "D"
17
18 timezone = None
19
20 fileIndex = None
21
22 blockIndex = None
23
24 path = None
25
26 #Hdf5 File
27
28 fpMetadata = None
29
30 listMetaname = None
31
32 listMetadata = None
33
34 fp = None
35
36 #dataOut reconstruction
37
38
39 dataOut = None
40
41 nChannels = None #Dimension 0
42
43 nPoints = None #Dimension 1, number of Points or Parameters
44
45 nSamples = None #Dimension 2, number of samples or ranges
46
47
48 def __init__(self):
49
50 return
51
52 def setup(self,path=None,
53 startDate=None,
54 endDate=None,
55 startTime=datetime.time(0,0,0),
56 endTime=datetime.time(23,59,59),
57 walk=True,
58 timezone='ut',
59 all=0,
60 online=False,
61 ext=None):
62
63 if ext==None:
64 ext = self.ext
65 self.timezone = timezone
66 # self.all = all
67 # self.online = online
68 self.path = path
69
70
71 if not(online):
72 #Busqueda de archivos offline
73 self.__searchFilesOffline(path, startDate, endDate, ext, startTime, endTime, walk)
74 else:
75 self.__searchFilesOnline(path, walk)
76
77 if not(self.filenameList):
78 print "There is no files into the folder: %s"%(path)
79 sys.exit(-1)
80
81 # self.__getExpParameters()
82
83 self.fileIndex = -1
84
85 self.__setNextFileOffline()
86
87 self.__readMetadata()
88
89 self.blockIndex = 0
90
91 return
92
93 def __searchFilesOffline(self,
94 path,
95 startDate,
96 endDate,
97 ext,
98 startTime=datetime.time(0,0,0),
99 endTime=datetime.time(23,59,59),
100 walk=True):
101
102 # self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
103 #
104 # self.__checkPath()
105 #
106 # self.__findDataForDates()
107 #
108 # self.__selectDataForTimes()
109 #
110 # for i in range(len(self.filenameList)):
111 # print "%s" %(self.filenameList[i])
112
113 pathList = []
114
115 if not walk:
116 #pathList.append(path)
117 multi_path = path.split(',')
118 for single_path in multi_path:
119 pathList.append(single_path)
120
121 else:
122 #dirList = []
123 multi_path = path.split(',')
124 for single_path in multi_path:
125 dirList = []
126 for thisPath in os.listdir(single_path):
127 if not os.path.isdir(os.path.join(single_path,thisPath)):
128 continue
129 if not isDoyFolder(thisPath):
130 continue
131
132 dirList.append(thisPath)
133
134 if not(dirList):
135 return None, None
136
137 thisDate = startDate
138
139 while(thisDate <= endDate):
140 year = thisDate.timetuple().tm_year
141 doy = thisDate.timetuple().tm_yday
142
143 matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*')
144 if len(matchlist) == 0:
145 thisDate += datetime.timedelta(1)
146 continue
147 for match in matchlist:
148 pathList.append(os.path.join(single_path,match))
149
150 thisDate += datetime.timedelta(1)
151
152 if pathList == []:
153 print "Any folder was found for the date range: %s-%s" %(startDate, endDate)
154 return None, None
155
156 print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate)
157
158 filenameList = []
159 datetimeList = []
160 pathDict = {}
161 filenameList_to_sort = []
162
163 for i in range(len(pathList)):
164
165 thisPath = pathList[i]
166
167 fileList = glob.glob1(thisPath, "*%s" %ext)
168 fileList.sort()
169 pathDict.setdefault(fileList[0])
170 pathDict[fileList[0]] = i
171 filenameList_to_sort.append(fileList[0])
172
173 filenameList_to_sort.sort()
174
175 for file in filenameList_to_sort:
176 thisPath = pathList[pathDict[file]]
177
178 fileList = glob.glob1(thisPath, "*%s" %ext)
179 fileList.sort()
180
181 for file in fileList:
182
183 filename = os.path.join(thisPath,file)
184 thisDatetime = self.__isFileinThisTime(filename, startTime, endTime)
185
186 if not(thisDatetime):
187 continue
188
189 filenameList.append(filename)
190 datetimeList.append(thisDatetime)
191
192 if not(filenameList):
193 print "Any file was found for the time range %s - %s" %(startTime, endTime)
194 return None, None
195
196 print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime)
197 print
198
199 for i in range(len(filenameList)):
200 print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
201
202 self.filenameList = filenameList
203 self.datetimeList = datetimeList
204
205 return pathList, filenameList
206
207 def __isFileinThisTime(self, filename, startTime, endTime):
208 """
209 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
210
211 Inputs:
212 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
213
214 startTime : tiempo inicial del rango seleccionado en formato datetime.time
215
216 endTime : tiempo final del rango seleccionado en formato datetime.time
217
218 Return:
219 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
220 fecha especificado, de lo contrario retorna False.
221
222 Excepciones:
223 Si el archivo no existe o no puede ser abierto
224 Si la cabecera no puede ser leida.
225
226 """
227
228
229 try:
230 fp = fp = h5py.File(filename,'r')
231 except IOError:
232 traceback.print_exc()
233 raise IOError, "The file %s can't be opened" %(filename)
234
235 grp = fp['Data']
236 time = grp['time']
237 time0 = time[:][0]
238
239 fp.close()
240
241 thisDatetime = datetime.datetime.utcfromtimestamp(time0)
242
243 if self.timezone == 'lt':
244 thisDatetime = thisDatetime - datetime.timedelta(minutes = 300)
245
246 thisTime = thisDatetime.time()
247
248 if not ((startTime <= thisTime) and (endTime > thisTime)):
249 return None
250
251 return thisDatetime
252
253 def __checkPath(self):
254 if os.path.exists(self.path):
255 self.status = 1
256 else:
257 self.status = 0
258 print 'Path:%s does not exists'%self.path
259
260 return
261
262 def __setNextFileOffline(self):
263 idFile = self.fileIndex
264 idFile += 1
265
266 if not(idFile < len(self.filenameList)):
267 self.flagNoMoreFiles = 1
268 print "No more Files"
269 return 0
270
271 filename = self.filenameList[idFile]
272
273 filePointer = h5py.File(filename,'r')
274
275 self.flagIsNewFile = 1
276 self.fileIndex = idFile
277 self.filename = filename
278
279 self.fp = filePointer
280
281 print "Setting the file: %s"%self.filename
282
283 self.__readMetadata()
284
285 return 1
286
287 def __readMetadata(self):
288 grp = self.fp['Data']
289 self.pathMeta = os.path.join(self.path, grp.attrs['metadata'])
290 filePointer = h5py.File(self.pathMeta,'r')
291 groupPointer = filePointer['Metadata']
292
293 listMetaname = []
294 listMetadata = []
295 for item in groupPointer.items():
296 name = item[0]
297
298 if name=='data shape':
299 self.nSamples = 1
300 self.nPoints = 1
301 self.nChannels = 1
302 else:
303 data = groupPointer[name][:]
304 listMetaname.append(name)
305 listMetadata.append(data)
306
307 if name=='type':
308 self.__initDataOut(name)
309
310 filePointer.close()
311
312 self.listMetadata = listMetaname
313 self.listMetadata = listMetadata
314
315 return
316
317 def __initDataOut(self, type):
318
319 if 'type'=='Parameters':
320 self.dataOut = Parameters()
321 elif 'type'=='Spectra':
322 self.dataOut = Spectra()
323 elif 'type'=='Voltage':
324 self.dataOut = Voltage()
325 elif 'type'=='Correlation':
326 self.dataOut = Correlation()
327
328 return
329
330 def __setDataOut(self):
331 listMetadata = self.listMetadata
332 listMetaname = self.listMetaname
333 listDataname = self.listDataname
334 listData = self.listData
335
336 blockIndex = self.blockIndex
337
338 for i in range(len(listMetadata)):
339 setattr(self.dataOut,listMetaname[i],listMetadata[i])
340
341 for j in range(len(listData)):
342 setattr(self.dataOut,listDataname[j][blockIndex,:],listData[j][blockIndex,:])
343
344 return
345
346 def getData(self):
347
348 if self.flagNoMoreFiles:
349 self.dataOut.flagNoData = True
350 print 'Process finished'
351 return 0
352
353 if self.__hasNotDataInBuffer():
354 self.__setNextFile()
355
356
357 if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
358 self.dataOut.flagNoData = True
359 return 0
360
361 self.__setDataOut()
362 self.dataOut.flagNoData = False
363
364 self.blockIndex += 1
365
366 return self.dataOut.data
367
368 def run(self, **kwargs):
369
370 if not(self.isConfig):
371 self.setup(**kwargs)
372 self.setObjProperties()
373 self.isConfig = True
374
375 self.getData()
376
377 return
378
379 class HDF5Writer(Operation):
380
381 ext = ".hdf5"
382
383 optchar = "D"
384
385 metaoptchar = "M"
386
387 metaFile = None
388
389 path = None
390
391 setFile = None
392
393 fp = None
394
395 grp = None
396
397 ds = None
398
399 firsttime = True
400
401 #Configurations
402
403 blocksPerFile = None
404
405 blockIndex = None
406
407 dataOut = None
408
409 #Data Arrays
410
411 dataList = None
412
413 metadataList = None
414
415 dataDim = None
416
417 def __init__(self):
418
419 Operation.__init__(self)
420 self.isConfig = False
421 return
422
423
424 def setup(self, dataOut, **kwargs):
425
426 self.path = kwargs['path']
427
428 if kwargs.has_key('ext'):
429 self.ext = kwargs['ext']
430 else:
431 self.blocksPerFile = 10
432
433 if kwargs.has_key('blocksPerFile'):
434 self.blocksPerFile = kwargs['blocksPerFile']
435 else:
436 self.blocksPerFile = 10
437
438 self.dataOut = dataOut
439
440 self.metadataList = ['inputUnit','abscissaRange','heightRange']
441
442 self.dataList = ['data_param', 'data_error', 'data_SNR']
443
444 self.dataDim = numpy.zeros((len(self.dataList),3))
445
446 for i in range(len(self.dataList)):
447
448 dataDim = getattr(self.dataOut, self.dataList[i]).shape
449
450 if len(dataDim) == 3:
451 self.dataDim[i,:] = numpy.array(dataDim)
452 else:
453 self.dataDim[i,:-1] = numpy.array(dataDim)
454 self.dataDim[i,-1] = numpy.nan
455
456 self.blockIndex = 0
457
458 return
459
460 def putMetadata(self):
461
462 fp = self.createMetadataFile()
463 self.writeMetadata(fp)
464 fp.close()
465 return
466
467 def createMetadataFile(self):
468 ext = self.ext
469 path = self.path
470 setFile = self.setFile
471
472 timeTuple = time.localtime(self.dataOut.utctime)
473 subfolder = ''
474
475 fullpath = os.path.join( path, subfolder )
476 if not( os.path.exists(fullpath) ):
477 os.mkdir(fullpath)
478 setFile = -1 #inicializo mi contador de seteo
479 else:
480 filesList = os.listdir( fullpath )
481 if len( filesList ) > 0:
482 filesList = sorted( filesList, key=str.lower )
483 filen = filesList[-1]
484 # el filename debera tener el siguiente formato
485 # 0 1234 567 89A BCDE (hex)
486 # x YYYY DDD SSS .ext
487 if isNumber( filen[8:11] ):
488 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
489 else:
490 setFile = -1
491 else:
492 setFile = -1 #inicializo mi contador de seteo
493
494 setFile += 1
495
496 file = '%s%4.4d%3.3d%3.3d%s' % (self.metaoptchar,
497 timeTuple.tm_year,
498 timeTuple.tm_yday,
499 setFile,
500 ext )
501
502 filename = os.path.join( path, subfolder, file )
503 self.metaFile = file
504 #Setting HDF5 File
505 fp = h5py.File(filename,'w')
506
507 return fp
508
509 def writeMetadata(self, fp):
510
511 grp = fp.create_group("Metadata")
512
513 for i in range(len(self.metadataList)):
514 grp.create_dataset(self.metadataList[i], data=getattr(self.dataOut, self.metadataList[i]))
515 return
516
517 def setNextFile(self):
518
519 ext = self.ext
520 path = self.path
521 setFile = self.setFile
522
523 if self.fp != None:
524 self.fp.close()
525
526 timeTuple = time.localtime(self.dataOut.utctime)
527 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
528
529 fullpath = os.path.join( path, subfolder )
530 if not( os.path.exists(fullpath) ):
531 os.mkdir(fullpath)
532 setFile = -1 #inicializo mi contador de seteo
533 else:
534 filesList = os.listdir( fullpath )
535 if len( filesList ) > 0:
536 filesList = sorted( filesList, key=str.lower )
537 filen = filesList[-1]
538 # el filename debera tener el siguiente formato
539 # 0 1234 567 89A BCDE (hex)
540 # x YYYY DDD SSS .ext
541 if isNumber( filen[8:11] ):
542 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
543 else:
544 setFile = -1
545 else:
546 setFile = -1 #inicializo mi contador de seteo
547
548 setFile += 1
549
550 file = '%s%4.4d%3.3d%3.3d%s' % (self.optchar,
551 timeTuple.tm_year,
552 timeTuple.tm_yday,
553 setFile,
554 ext )
555
556 filename = os.path.join( path, subfolder, file )
557
558 #Setting HDF5 File
559 fp = h5py.File(filename,'w')
560 grp = fp.create_group("Data")
561 grp.attrs['metadata'] = self.metaFile
562
563
564
565 ds = []
566 data = []
567
568 for i in range(len(self.dataList)):
569
570 grp0 = grp.create_group(self.dataList[i])
571
572 for j in range(int(self.dataDim[i,0])):
573 tableName = "channel" + str(j)
574
575 if not(numpy.isnan(self.dataDim[i,2])):
576 ds0 = grp0.create_dataset(tableName, (1,1,1) , chunks = True)
577 else:
578 ds0 = grp0.create_dataset(tableName, (1,1) , chunks = True)
579
580 ds.append(ds0)
581 data.append([])
582
583 ds0 = grp.create_dataset("time", (1,) , chunks = True)
584 ds.append(ds0)
585 data.append([])
586
587 #Saving variables
588 print 'Writing the file: %s'%filename
589 self.fp = fp
590 self.grp = grp
591 self.ds = ds
592 self.data = data
593
594 self.setFile = setFile
595 self.firsttime = True
596 self.blockIndex = 0
597 return
598
599 def putData(self):
600 self.setBlock()
601 self.writeBlock()
602
603 if self.blockIndex == self.blocksPerFile:
604 self.setNextFile()
605 return
606
607 def setBlock(self):
608
609 #Creating Arrays
610 data = self.data
611 ind = 0
612 for i in range(len(self.dataList)):
613 dataAux = getattr(self.dataOut,self.dataList[i])
614
615 for j in range(int(self.dataDim[i,0])):
616 data[ind] = dataAux[j,:]
617 if not(numpy.isnan(self.dataDim[i,2])):
618 data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
619 if not self.firsttime:
620 data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
621 else:
622 data[ind] = data[ind].reshape((1,data[ind].shape[0]))
623 if not self.firsttime:
624 data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
625 ind += 1
626
627 data[ind] = numpy.array([self.dataOut.utctime])
628 if not self.firsttime:
629 self.data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
630 self.data = data
631
632 return
633
634 def writeBlock(self):
635
636 for i in range(len(self.ds)):
637 self.ds[i].shape = self.data[i].shape
638 self.ds[i][:] = self.data[i]
639
640 self.blockIndex += 1
641
642 self.grp['blocksPerFile'] = self.blockIndex
643
644 self.firsttime = False
645 return
646
647 def run(self, dataOut, **kwargs):
648 if not(self.isConfig):
649 self.setup(dataOut, **kwargs)
650 self.isConfig = True
651 self.putMetadata()
652 self.setNextFile()
653
654 self.putData()
655 return
656
Show file before | Show file after | Hide comments
@@ -0,0 +1,65
1 # DIAS 19 Y 20 FEB 2014
2 # Comprobacion de Resultados DBS con SA
3
4 import os, sys
5
6 path = os.path.split(os.getcwd())[0]
7 sys.path.append(path)
8
9 from controller import *
10
11 desc = "DBS Experiment Test"
12 filename = "DBStest.xml"
13
14 controllerObj = Project()
15
16 controllerObj.setup(id = '191', name='test01', description=desc)
17
18 #Experimentos
19
20 path = "/home/soporte/Data/drifts/HDF5"
21 pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba'
22 pathFile = '/home/soporte/Data/drifts/HDF5'
23
24 xmin = 0
25 xmax = 24
26 #------------------------------------------------------------------------------------------------
27 readUnitConfObj = controllerObj.addReadUnit(datatype='HDF5Reader',
28 path=path,
29 startDate='2012/09/06',
30 endDate='2012/09/06',
31 startTime='00:00:00',
32 endTime='23:59:59',
33 timezone='lt',
34 walk=1)
35
36 procUnitConfObj0 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=readUnitConfObj.getId())
37 #--------------------------------------------------------------------------------------------------
38
39 opObj11 = procUnitConfObj0.addOperation(name='EWDriftsEstimation', optype='other')
40 opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist')
41 opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float')
42
43 opObj23 = procUnitConfObj0.addOperation(name='EWDriftsPlot', optype='other')
44 opObj23.addParameter(name='id', value='4', format='int')
45 opObj23.addParameter(name='wintitle', value='EW Drifts', format='str')
46 opObj23.addParameter(name='save', value='1', format='bool')
47 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
48 opObj23.addParameter(name='zminZonal', value='-150', format='int')
49 opObj23.addParameter(name='zmaxZonal', value='150', format='int')
50 opObj23.addParameter(name='zminVertical', value='-30', format='float')
51 opObj23.addParameter(name='zmaxVertical', value='30', format='float')
52 opObj23.addParameter(name='SNR_1', value='1', format='bool')
53 opObj23.addParameter(name='SNRmax', value='5', format='int')
54 # opObj23.addParameter(name='SNRthresh', value='-50', format='float')
55 opObj23.addParameter(name='xmin', value=xmin, format='float')
56 opObj23.addParameter(name='xmax', value=xmax, format='float')
57 #--------------------------------------------------------------------------------------------------
58 print "Escribiendo el archivo XML"
59 controllerObj.writeXml(filename)
60 print "Leyendo el archivo XML"
61 controllerObj.readXml(filename)
62
63 controllerObj.createObjects()
64 controllerObj.connectObjects()
65 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,746 +1,756
1 1 '''
2 2 Created on September , 2012
3 3 @author:
4 4 '''
5 5 from xml.etree.ElementTree import Element, SubElement, ElementTree
6 6 from xml.etree import ElementTree as ET
7 7 from xml.dom import minidom
8 8
9 9 import datetime
10 10 from model import *
11 11
12 import ast
13
12 14 def prettify(elem):
13 15 """Return a pretty-printed XML string for the Element.
14 16 """
15 17 rough_string = ET.tostring(elem, 'utf-8')
16 18 reparsed = minidom.parseString(rough_string)
17 19 return reparsed.toprettyxml(indent=" ")
18 20
19 21 class ParameterConf():
20 22
21 23 id = None
22 24 name = None
23 25 value = None
24 26 format = None
25 27
26 28 __value = None
27 29
28 30 ELEMENTNAME = 'Parameter'
29 31
30 32 def __init__(self):
31 33
32 34 self.format = 'str'
33 35
34 36 def getElementName(self):
35 37
36 38 return self.ELEMENTNAME
37 39
38 40 def getValue(self):
39 41
40 42 if self.__value != None:
41 43
42 44 return self.__value
43 45
44 46 value = self.value
45 47
46 48 if self.format == 'list':
47 49 strList = value.split(',')
48 50 return strList
49 51
50 52 if self.format == 'intlist':
51 53 strList = value.split(',')
52 54 intList = [int(x) for x in strList]
53 55 return intList
54 56
55 57 if self.format == 'floatlist':
56 58 strList = value.split(',')
57 59 floatList = [float(x) for x in strList]
58 60 return floatList
59 61
60 62 if self.format == 'date':
61 63 strList = value.split('/')
62 64 intList = [int(x) for x in strList]
63 65 date = datetime.date(intList[0], intList[1], intList[2])
64 66 return date
65 67
66 68 if self.format == 'time':
67 69 strList = value.split(':')
68 70 intList = [int(x) for x in strList]
69 71 time = datetime.time(intList[0], intList[1], intList[2])
70 72 return time
71 73
72 74 if self.format == 'bool':
73 75 value = int(value)
74 76
75 77 if self.format == 'pairsList':
76 78 """
77 79 Example:
78 80 value = (0,1),(1,2)
79 81 """
80 82
81 83 value = value.replace('(', '')
82 84 value = value.replace(')', '')
83 85
84 86 strList = value.split(',')
85 87 intList = [int(item) for item in strList]
86 88 pairList = []
87 89 for i in range(len(intList)/2):
88 90 pairList.append((intList[i*2], intList[i*2 + 1]))
89 91
90 92 return pairList
91 93
94 if self.format == 'multiList':
95 """
96 Example:
97 value = (0,1,2),(3,4,5)
98 """
99 multiList = ast.literal_eval(value)
100 return multiList
101
92 102 func = eval(self.format)
93 103
94 104 self.__value = func(value)
95 105
96 106 return self.__value
97 107
98 108 def setup(self, id, name, value, format='str'):
99 109
100 110 self.id = id
101 111 self.name = name
102 112 self.value = str(value)
103 113 self.format = format
104 114
105 115 def makeXml(self, opElement):
106 116
107 117 parmElement = SubElement(opElement, self.ELEMENTNAME)
108 118 parmElement.set('id', str(self.id))
109 119 parmElement.set('name', self.name)
110 120 parmElement.set('value', self.value)
111 121 parmElement.set('format', self.format)
112 122
113 123 def readXml(self, parmElement):
114 124
115 125 self.id = parmElement.get('id')
116 126 self.name = parmElement.get('name')
117 127 self.value = parmElement.get('value')
118 128 self.format = parmElement.get('format')
119 129
120 130 def printattr(self):
121 131
122 132 print "Parameter[%s]: name = %s, value = %s, format = %s" %(self.id, self.name, self.value, self.format)
123 133
124 134 class OperationConf():
125 135
126 136 id = None
127 137 name = None
128 138 priority = None
129 139 type = None
130 140
131 141 parmConfObjList = []
132 142
133 143 ELEMENTNAME = 'Operation'
134 144
135 145 def __init__(self):
136 146
137 147 id = 0
138 148 name = None
139 149 priority = None
140 150 type = 'self'
141 151
142 152
143 153 def __getNewId(self):
144 154
145 155 return int(self.id)*10 + len(self.parmConfObjList) + 1
146 156
147 157 def getElementName(self):
148 158
149 159 return self.ELEMENTNAME
150 160
151 161 def getParameterObjList(self):
152 162
153 163 return self.parmConfObjList
154 164
155 165 def setup(self, id, name, priority, type):
156 166
157 167 self.id = id
158 168 self.name = name
159 169 self.type = type
160 170 self.priority = priority
161 171
162 172 self.parmConfObjList = []
163 173
164 174 def addParameter(self, name, value, format='str'):
165 175
166 176 id = self.__getNewId()
167 177
168 178 parmConfObj = ParameterConf()
169 179 parmConfObj.setup(id, name, value, format)
170 180
171 181 self.parmConfObjList.append(parmConfObj)
172 182
173 183 return parmConfObj
174 184
175 185 def makeXml(self, upElement):
176 186
177 187 opElement = SubElement(upElement, self.ELEMENTNAME)
178 188 opElement.set('id', str(self.id))
179 189 opElement.set('name', self.name)
180 190 opElement.set('type', self.type)
181 191 opElement.set('priority', str(self.priority))
182 192
183 193 for parmConfObj in self.parmConfObjList:
184 194 parmConfObj.makeXml(opElement)
185 195
186 196 def readXml(self, opElement):
187 197
188 198 self.id = opElement.get('id')
189 199 self.name = opElement.get('name')
190 200 self.type = opElement.get('type')
191 201 self.priority = opElement.get('priority')
192 202
193 203 self.parmConfObjList = []
194 204
195 205 parmElementList = opElement.getiterator(ParameterConf().getElementName())
196 206
197 207 for parmElement in parmElementList:
198 208 parmConfObj = ParameterConf()
199 209 parmConfObj.readXml(parmElement)
200 210 self.parmConfObjList.append(parmConfObj)
201 211
202 212 def printattr(self):
203 213
204 214 print "%s[%s]: name = %s, type = %s, priority = %s" %(self.ELEMENTNAME,
205 215 self.id,
206 216 self.name,
207 217 self.type,
208 218 self.priority)
209 219
210 220 for parmConfObj in self.parmConfObjList:
211 221 parmConfObj.printattr()
212 222
213 223 def createObject(self):
214 224
215 225 if self.type == 'self':
216 226 raise ValueError, "This operation type cannot be created"
217 227
218 228 if self.type == 'external' or self.type == 'other':
219 229 className = eval(self.name)
220 230 opObj = className()
221 231
222 232 return opObj
223 233
224 234 class ProcUnitConf():
225 235
226 236 id = None
227 237 name = None
228 238 datatype = None
229 239 inputId = None
230 240
231 241 opConfObjList = []
232 242
233 243 procUnitObj = None
234 244 opObjList = []
235 245
236 246 ELEMENTNAME = 'ProcUnit'
237 247
238 248 def __init__(self):
239 249
240 250 self.id = None
241 251 self.datatype = None
242 252 self.name = None
243 253 self.inputId = None
244 254
245 255 self.opConfObjList = []
246 256
247 257 self.procUnitObj = None
248 258 self.opObjDict = {}
249 259
250 260 def __getPriority(self):
251 261
252 262 return len(self.opConfObjList)+1
253 263
254 264 def __getNewId(self):
255 265
256 266 return int(self.id)*10 + len(self.opConfObjList) + 1
257 267
258 268 def getElementName(self):
259 269
260 270 return self.ELEMENTNAME
261 271
262 272 def getId(self):
263 273
264 274 return str(self.id)
265 275
266 276 def getInputId(self):
267 277
268 278 return str(self.inputId)
269 279
270 280 def getOperationObjList(self):
271 281
272 282 return self.opConfObjList
273 283
274 284 def getProcUnitObj(self):
275 285
276 286 return self.procUnitObj
277 287
278 288 def setup(self, id, name, datatype, inputId):
279 289
280 290 self.id = id
281 291 self.name = name
282 292 self.datatype = datatype
283 293 self.inputId = inputId
284 294
285 295 self.opConfObjList = []
286 296
287 297 self.addOperation(name='run', optype='self')
288 298
289 299 def addParameter(self, **kwargs):
290 300
291 301 opObj = self.opConfObjList[0]
292 302
293 303 opObj.addParameter(**kwargs)
294 304
295 305 return opObj
296 306
297 307 def addOperation(self, name, optype='self'):
298 308
299 309 id = self.__getNewId()
300 310 priority = self.__getPriority()
301 311
302 312 opConfObj = OperationConf()
303 313 opConfObj.setup(id, name=name, priority=priority, type=optype)
304 314
305 315 self.opConfObjList.append(opConfObj)
306 316
307 317 return opConfObj
308 318
309 319 def makeXml(self, procUnitElement):
310 320
311 321 upElement = SubElement(procUnitElement, self.ELEMENTNAME)
312 322 upElement.set('id', str(self.id))
313 323 upElement.set('name', self.name)
314 324 upElement.set('datatype', self.datatype)
315 325 upElement.set('inputId', str(self.inputId))
316 326
317 327 for opConfObj in self.opConfObjList:
318 328 opConfObj.makeXml(upElement)
319 329
320 330 def readXml(self, upElement):
321 331
322 332 self.id = upElement.get('id')
323 333 self.name = upElement.get('name')
324 334 self.datatype = upElement.get('datatype')
325 335 self.inputId = upElement.get('inputId')
326 336
327 337 self.opConfObjList = []
328 338
329 339 opElementList = upElement.getiterator(OperationConf().getElementName())
330 340
331 341 for opElement in opElementList:
332 342 opConfObj = OperationConf()
333 343 opConfObj.readXml(opElement)
334 344 self.opConfObjList.append(opConfObj)
335 345
336 346 def printattr(self):
337 347
338 348 print "%s[%s]: name = %s, datatype = %s, inputId = %s" %(self.ELEMENTNAME,
339 349 self.id,
340 350 self.name,
341 351 self.datatype,
342 352 self.inputId)
343 353
344 354 for opConfObj in self.opConfObjList:
345 355 opConfObj.printattr()
346 356
347 357 def createObjects(self):
348 358
349 359 className = eval(self.name)
350 360 procUnitObj = className()
351 361
352 362 for opConfObj in self.opConfObjList:
353 363
354 364 if opConfObj.type == 'self':
355 365 continue
356 366
357 367 opObj = opConfObj.createObject()
358 368
359 369 self.opObjDict[opConfObj.id] = opObj
360 370 procUnitObj.addOperation(opObj, opConfObj.id)
361 371
362 372 self.procUnitObj = procUnitObj
363 373
364 374 return procUnitObj
365 375
366 376 def run(self):
367 377
368 378 finalSts = False
369 379
370 380 for opConfObj in self.opConfObjList:
371 381
372 382 kwargs = {}
373 383 for parmConfObj in opConfObj.getParameterObjList():
374 384 kwargs[parmConfObj.name] = parmConfObj.getValue()
375 385
376 386 #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id)
377 387 sts = self.procUnitObj.call(opType = opConfObj.type,
378 388 opName = opConfObj.name,
379 389 opId = opConfObj.id,
380 390 **kwargs)
381 391 finalSts = finalSts or sts
382 392
383 393 return finalSts
384 394
385 395 class ReadUnitConf(ProcUnitConf):
386 396
387 397 path = None
388 398 startDate = None
389 399 endDate = None
390 400 startTime = None
391 401 endTime = None
392 402
393 403 ELEMENTNAME = 'ReadUnit'
394 404
395 405 def __init__(self):
396 406
397 407 self.id = None
398 408 self.datatype = None
399 409 self.name = None
400 410 self.inputId = 0
401 411
402 412 self.opConfObjList = []
403 413 self.opObjList = []
404 414
405 415 def getElementName(self):
406 416
407 417 return self.ELEMENTNAME
408 418
409 419 def setup(self, id, name, datatype, path="", startDate="", endDate="", startTime="", endTime="", **kwargs):
410 420
411 421 self.id = id
412 422 self.name = name
413 423 self.datatype = datatype
414 424
415 425 self.path = path
416 426 self.startDate = startDate
417 427 self.endDate = endDate
418 428 self.startTime = startTime
419 429 self.endTime = endTime
420 430
421 431 self.addRunOperation(**kwargs)
422 432
423 433 def addRunOperation(self, **kwargs):
424 434
425 435 opObj = self.addOperation(name = 'run', optype = 'self')
426 436
427 437 opObj.addParameter(name='path' , value=self.path, format='str')
428 438 opObj.addParameter(name='startDate' , value=self.startDate, format='date')
429 439 opObj.addParameter(name='endDate' , value=self.endDate, format='date')
430 440 opObj.addParameter(name='startTime' , value=self.startTime, format='time')
431 441 opObj.addParameter(name='endTime' , value=self.endTime, format='time')
432 442
433 443 for key, value in kwargs.items():
434 444 opObj.addParameter(name=key, value=value, format=type(value).__name__)
435 445
436 446 return opObj
437 447
438 448
439 449 class Project():
440 450
441 451 id = None
442 452 name = None
443 453 description = None
444 454 # readUnitConfObjList = None
445 455 procUnitConfObjDict = None
446 456
447 457 ELEMENTNAME = 'Project'
448 458
449 459 def __init__(self):
450 460
451 461 self.id = None
452 462 self.name = None
453 463 self.description = None
454 464
455 465 # self.readUnitConfObjList = []
456 466 self.procUnitConfObjDict = {}
457 467
458 468 def __getNewId(self):
459 469
460 470 id = int(self.id)*10 + len(self.procUnitConfObjDict) + 1
461 471
462 472 return str(id)
463 473
464 474 def getElementName(self):
465 475
466 476 return self.ELEMENTNAME
467 477
468 478 def setup(self, id, name, description):
469 479
470 480 self.id = id
471 481 self.name = name
472 482 self.description = description
473 483
474 484 def addReadUnit(self, datatype, **kwargs):
475 485
476 486 id = self.__getNewId()
477 487 name = '%s' %(datatype)
478 488
479 489 readUnitConfObj = ReadUnitConf()
480 490 readUnitConfObj.setup(id, name, datatype, **kwargs)
481 491
482 492 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
483 493
484 494 return readUnitConfObj
485 495
486 496 def addProcUnit(self, datatype, inputId):
487 497
488 498 id = self.__getNewId()
489 499 name = '%s' %(datatype)
490 500
491 501 procUnitConfObj = ProcUnitConf()
492 502 procUnitConfObj.setup(id, name, datatype, inputId)
493 503
494 504 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
495 505
496 506 return procUnitConfObj
497 507
498 508 def makeXml(self):
499 509
500 510 projectElement = Element('Project')
501 511 projectElement.set('id', str(self.id))
502 512 projectElement.set('name', self.name)
503 513 projectElement.set('description', self.description)
504 514
505 515 # for readUnitConfObj in self.readUnitConfObjList:
506 516 # readUnitConfObj.makeXml(projectElement)
507 517
508 518 for procUnitConfObj in self.procUnitConfObjDict.values():
509 519 procUnitConfObj.makeXml(projectElement)
510 520
511 521 self.projectElement = projectElement
512 522
513 523 def writeXml(self, filename):
514 524
515 525 self.makeXml()
516 526
517 527 print prettify(self.projectElement)
518 528
519 529 ElementTree(self.projectElement).write(filename, method='xml')
520 530
521 531 def readXml(self, filename):
522 532
523 533 #tree = ET.parse(filename)
524 534 self.projectElement = None
525 535 # self.readUnitConfObjList = []
526 536 self.procUnitConfObjDict = {}
527 537
528 538 self.projectElement = ElementTree().parse(filename)
529 539
530 540 self.project = self.projectElement.tag
531 541
532 542 self.id = self.projectElement.get('id')
533 543 self.name = self.projectElement.get('name')
534 544 self.description = self.projectElement.get('description')
535 545
536 546 readUnitElementList = self.projectElement.getiterator(ReadUnitConf().getElementName())
537 547
538 548 for readUnitElement in readUnitElementList:
539 549 readUnitConfObj = ReadUnitConf()
540 550 readUnitConfObj.readXml(readUnitElement)
541 551
542 552 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
543 553
544 554 procUnitElementList = self.projectElement.getiterator(ProcUnitConf().getElementName())
545 555
546 556 for procUnitElement in procUnitElementList:
547 557 procUnitConfObj = ProcUnitConf()
548 558 procUnitConfObj.readXml(procUnitElement)
549 559
550 560 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
551 561
552 562 def printattr(self):
553 563
554 564 print "Project[%s]: name = %s, description = %s" %(self.id,
555 565 self.name,
556 566 self.description)
557 567
558 568 # for readUnitConfObj in self.readUnitConfObjList:
559 569 # readUnitConfObj.printattr()
560 570
561 571 for procUnitConfObj in self.procUnitConfObjDict.values():
562 572 procUnitConfObj.printattr()
563 573
564 574 def createObjects(self):
565 575
566 576 # for readUnitConfObj in self.readUnitConfObjList:
567 577 # readUnitConfObj.createObjects()
568 578
569 579 for procUnitConfObj in self.procUnitConfObjDict.values():
570 580 procUnitConfObj.createObjects()
571 581
572 582 def __connect(self, objIN, thisObj):
573 583
574 584 thisObj.setInput(objIN.getOutputObj())
575 585
576 586 def connectObjects(self):
577 587
578 588 for thisPUConfObj in self.procUnitConfObjDict.values():
579 589
580 590 inputId = thisPUConfObj.getInputId()
581 591
582 592 if int(inputId) == 0:
583 593 continue
584 594
585 595 #Get input object
586 596 puConfINObj = self.procUnitConfObjDict[inputId]
587 597 puObjIN = puConfINObj.getProcUnitObj()
588 598
589 599 #Get current object
590 600 thisPUObj = thisPUConfObj.getProcUnitObj()
591 601
592 602 self.__connect(puObjIN, thisPUObj)
593 603
594 604 def run(self):
595 605
596 606 # for readUnitConfObj in self.readUnitConfObjList:
597 607 # readUnitConfObj.run()
598 608
599 609 while(True):
600 610
601 611 finalSts = False
602 612
603 613 for procUnitConfObj in self.procUnitConfObjDict.values():
604 614 #print "Running the '%s' process with %s" %(procUnitConfObj.name, procUnitConfObj.id)
605 615 sts = procUnitConfObj.run()
606 616 finalSts = finalSts or sts
607 617
608 618 #If every process unit finished so end process
609 619 if not(finalSts):
610 620 print "Every process units have finished"
611 621 break
612 622
613 623 if __name__ == '__main__':
614 624
615 625 desc = "Segundo Test"
616 626 filename = "schain.xml"
617 627
618 628 controllerObj = Project()
619 629
620 630 controllerObj.setup(id = '191', name='test01', description=desc)
621 631
622 632 readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage',
623 633 path='data/rawdata/',
624 634 startDate='2011/01/01',
625 635 endDate='2012/12/31',
626 636 startTime='00:00:00',
627 637 endTime='23:59:59',
628 638 online=1,
629 639 walk=1)
630 640
631 641 # opObj00 = readUnitConfObj.addOperation(name='printInfo')
632 642
633 643 procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId())
634 644
635 645 opObj10 = procUnitConfObj0.addOperation(name='selectChannels')
636 646 opObj10.addParameter(name='channelList', value='3,4,5', format='intlist')
637 647
638 648 opObj10 = procUnitConfObj0.addOperation(name='selectHeights')
639 649 opObj10.addParameter(name='minHei', value='90', format='float')
640 650 opObj10.addParameter(name='maxHei', value='180', format='float')
641 651
642 652 opObj12 = procUnitConfObj0.addOperation(name='CohInt', optype='external')
643 653 opObj12.addParameter(name='n', value='10', format='int')
644 654
645 655 procUnitConfObj1 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId())
646 656 procUnitConfObj1.addParameter(name='nFFTPoints', value='32', format='int')
647 657 # procUnitConfObj1.addParameter(name='pairList', value='(0,1),(0,2),(1,2)', format='')
648 658
649 659
650 660 opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external')
651 661 opObj11.addParameter(name='idfigure', value='1', format='int')
652 662 opObj11.addParameter(name='wintitle', value='SpectraPlot0', format='str')
653 663 opObj11.addParameter(name='zmin', value='40', format='int')
654 664 opObj11.addParameter(name='zmax', value='90', format='int')
655 665 opObj11.addParameter(name='showprofile', value='1', format='int')
656 666
657 667 # opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='external')
658 668 # opObj11.addParameter(name='idfigure', value='2', format='int')
659 669 # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot', format='str')
660 670 # opObj11.addParameter(name='zmin', value='40', format='int')
661 671 # opObj11.addParameter(name='zmax', value='90', format='int')
662 672
663 673
664 674 # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Voltage', inputId=procUnitConfObj0.getId())
665 675 #
666 676 # opObj12 = procUnitConfObj2.addOperation(name='CohInt', optype='external')
667 677 # opObj12.addParameter(name='n', value='2', format='int')
668 678 # opObj12.addParameter(name='overlapping', value='1', format='int')
669 679 #
670 680 # procUnitConfObj3 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj2.getId())
671 681 # procUnitConfObj3.addParameter(name='nFFTPoints', value='32', format='int')
672 682 #
673 683 # opObj11 = procUnitConfObj3.addOperation(name='SpectraPlot', optype='external')
674 684 # opObj11.addParameter(name='idfigure', value='2', format='int')
675 685 # opObj11.addParameter(name='wintitle', value='SpectraPlot1', format='str')
676 686 # opObj11.addParameter(name='zmin', value='40', format='int')
677 687 # opObj11.addParameter(name='zmax', value='90', format='int')
678 688 # opObj11.addParameter(name='showprofile', value='1', format='int')
679 689
680 690 # opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='external')
681 691 # opObj11.addParameter(name='idfigure', value='10', format='int')
682 692 # opObj11.addParameter(name='wintitle', value='RTI', format='str')
683 693 ## opObj11.addParameter(name='xmin', value='21', format='float')
684 694 ## opObj11.addParameter(name='xmax', value='22', format='float')
685 695 # opObj11.addParameter(name='zmin', value='40', format='int')
686 696 # opObj11.addParameter(name='zmax', value='90', format='int')
687 697 # opObj11.addParameter(name='showprofile', value='1', format='int')
688 698 # opObj11.addParameter(name='timerange', value=str(60), format='int')
689 699
690 700 # opObj10 = procUnitConfObj1.addOperation(name='selectChannels')
691 701 # opObj10.addParameter(name='channelList', value='0,2,4,6', format='intlist')
692 702 #
693 703 # opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='external')
694 704 # opObj12.addParameter(name='n', value='2', format='int')
695 705 #
696 706 # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external')
697 707 # opObj11.addParameter(name='idfigure', value='2', format='int')
698 708 # opObj11.addParameter(name='wintitle', value='SpectraPlot10', format='str')
699 709 # opObj11.addParameter(name='zmin', value='70', format='int')
700 710 # opObj11.addParameter(name='zmax', value='90', format='int')
701 711 #
702 712 # opObj10 = procUnitConfObj1.addOperation(name='selectChannels')
703 713 # opObj10.addParameter(name='channelList', value='2,6', format='intlist')
704 714 #
705 715 # opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='external')
706 716 # opObj12.addParameter(name='n', value='2', format='int')
707 717 #
708 718 # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external')
709 719 # opObj11.addParameter(name='idfigure', value='3', format='int')
710 720 # opObj11.addParameter(name='wintitle', value='SpectraPlot10', format='str')
711 721 # opObj11.addParameter(name='zmin', value='70', format='int')
712 722 # opObj11.addParameter(name='zmax', value='90', format='int')
713 723
714 724
715 725 # opObj12 = procUnitConfObj1.addOperation(name='decoder')
716 726 # opObj12.addParameter(name='ncode', value='2', format='int')
717 727 # opObj12.addParameter(name='nbauds', value='8', format='int')
718 728 # opObj12.addParameter(name='code0', value='001110011', format='int')
719 729 # opObj12.addParameter(name='code1', value='001110011', format='int')
720 730
721 731
722 732
723 733 # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj1.getId())
724 734 #
725 735 # opObj21 = procUnitConfObj2.addOperation(name='IncohInt', optype='external')
726 736 # opObj21.addParameter(name='n', value='2', format='int')
727 737 #
728 738 # opObj11 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='external')
729 739 # opObj11.addParameter(name='idfigure', value='4', format='int')
730 740 # opObj11.addParameter(name='wintitle', value='SpectraPlot OBJ 2', format='str')
731 741 # opObj11.addParameter(name='zmin', value='70', format='int')
732 742 # opObj11.addParameter(name='zmax', value='90', format='int')
733 743
734 744 print "Escribiendo el archivo XML"
735 745
736 746 controllerObj.writeXml(filename)
737 747
738 748 print "Leyendo el archivo XML"
739 749 controllerObj.readXml(filename)
740 750 #controllerObj.printattr()
741 751
742 752 controllerObj.createObjects()
743 753 controllerObj.connectObjects()
744 754 controllerObj.run()
745 755
746 756 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,983 +1,983
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52
53 53 data = data.copy()
54 54
55 55 sortdata = numpy.sort(data,axis=None)
56 56 lenOfData = len(sortdata)
57 57 nums_min = lenOfData/10
58 58
59 59 if (lenOfData/10) > 2:
60 60 nums_min = lenOfData/10
61 61 else:
62 62 nums_min = 2
63 63
64 64 sump = 0.
65 65
66 66 sumq = 0.
67 67
68 68 j = 0
69 69
70 70 cont = 1
71 71
72 72 while((cont==1)and(j<lenOfData)):
73 73
74 74 sump += sortdata[j]
75 75
76 76 sumq += sortdata[j]**2
77 77
78 78 j += 1
79 79
80 80 if j > nums_min:
81 81 rtest = float(j)/(j-1) + 1.0/navg
82 82 if ((sumq*j) > (rtest*sump**2)):
83 83 j = j - 1
84 84 sump = sump - sortdata[j]
85 85 sumq = sumq - sortdata[j]**2
86 86 cont = 0
87 87
88 88 lnoise = sump /j
89 89 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
90 90 return lnoise
91 91
92 92 class Beam:
93 93 def __init__(self):
94 94 self.codeList = []
95 95 self.azimuthList = []
96 96 self.zenithList = []
97 97
98 98 class GenericData(object):
99 99
100 100 flagNoData = True
101 101
102 102 def __init__(self):
103 103
104 104 raise ValueError, "This class has not been implemented"
105 105
106 106 def copy(self, inputObj=None):
107 107
108 108 if inputObj == None:
109 109 return copy.deepcopy(self)
110 110
111 111 for key in inputObj.__dict__.keys():
112 112 self.__dict__[key] = inputObj.__dict__[key]
113 113
114 114 def deepcopy(self):
115 115
116 116 return copy.deepcopy(self)
117 117
118 118 def isEmpty(self):
119 119
120 120 return self.flagNoData
121 121
122 122 class JROData(GenericData):
123 123
124 124 # m_BasicHeader = BasicHeader()
125 125 # m_ProcessingHeader = ProcessingHeader()
126 126
127 127 systemHeaderObj = SystemHeader()
128 128
129 129 radarControllerHeaderObj = RadarControllerHeader()
130 130
131 131 # data = None
132 132
133 133 type = None
134 134
135 135 datatype = None #dtype but in string
136 136
137 137 # dtype = None
138 138
139 139 # nChannels = None
140 140
141 141 # nHeights = None
142 142
143 143 nProfiles = None
144 144
145 145 heightList = None
146 146
147 147 channelList = None
148 148
149 149 flagTimeBlock = False
150 150
151 151 useLocalTime = False
152 152
153 153 utctime = None
154 154
155 155 timeZone = None
156 156
157 157 dstFlag = None
158 158
159 159 errorCount = None
160 160
161 161 blocksize = None
162 162
163 163 nCode = None
164 164
165 165 nBaud = None
166 166
167 167 code = None
168 168
169 169 flagDecodeData = False #asumo q la data no esta decodificada
170 170
171 171 flagDeflipData = False #asumo q la data no esta sin flip
172 172
173 173 flagShiftFFT = False
174 174
175 175 # ippSeconds = None
176 176
177 177 timeInterval = None
178 178
179 179 nCohInt = None
180 180
181 181 noise = None
182 182
183 183 windowOfFilter = 1
184 184
185 185 #Speed of ligth
186 186 C = 3e8
187 187
188 188 frequency = 49.92e6
189 189
190 190 realtime = False
191 191
192 192 beacon_heiIndexList = None
193 193
194 194 last_block = None
195 195
196 196 blocknow = None
197 197
198 198 azimuth = None
199 199
200 200 zenith = None
201 201
202 202 beam = Beam()
203 203
204 204 def __init__(self):
205 205
206 206 raise ValueError, "This class has not been implemented"
207 207
208 208 def getNoise(self):
209 209
210 210 raise ValueError, "Not implemented"
211 211
212 212 def getNChannels(self):
213 213
214 214 return len(self.channelList)
215 215
216 216 def getChannelIndexList(self):
217 217
218 218 return range(self.nChannels)
219 219
220 220 def getNHeights(self):
221 221
222 222 return len(self.heightList)
223 223
224 224 def getHeiRange(self, extrapoints=0):
225 225
226 226 heis = self.heightList
227 227 # deltah = self.heightList[1] - self.heightList[0]
228 228 #
229 229 # heis.append(self.heightList[-1])
230 230
231 231 return heis
232 232
233 233 def getltctime(self):
234 234
235 235 if self.useLocalTime:
236 236 return self.utctime - self.timeZone*60
237 237
238 238 return self.utctime
239 239
240 240 def getDatatime(self):
241 241
242 242 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
243 243 return datatimeValue
244 244
245 245 def getTimeRange(self):
246 246
247 247 datatime = []
248 248
249 249 datatime.append(self.ltctime)
250 250 datatime.append(self.ltctime + self.timeInterval)
251 251
252 252 datatime = numpy.array(datatime)
253 253
254 254 return datatime
255 255
256 256 def getFmax(self):
257 257
258 258 PRF = 1./(self.ippSeconds * self.nCohInt)
259 259
260 260 fmax = PRF/2.
261 261
262 262 return fmax
263 263
264 264 def getVmax(self):
265 265
266 266 _lambda = self.C/self.frequency
267 267
268 268 vmax = self.getFmax() * _lambda
269 269
270 270 return vmax
271 271
272 272 def get_ippSeconds(self):
273 273 '''
274 274 '''
275 275 return self.radarControllerHeaderObj.ippSeconds
276 276
277 277 def set_ippSeconds(self, ippSeconds):
278 278 '''
279 279 '''
280 280
281 281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282 282
283 283 return
284 284
285 285 def get_dtype(self):
286 286 '''
287 287 '''
288 288 return getNumpyDtype(self.datatype)
289 289
290 290 def set_dtype(self, numpyDtype):
291 291 '''
292 292 '''
293 293
294 294 self.datatype = getDataTypeCode(numpyDtype)
295 295
296 296 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
297 297 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
298 298 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
299 299 #noise = property(getNoise, "I'm the 'nHeights' property.")
300 300 datatime = property(getDatatime, "I'm the 'datatime' property")
301 301 ltctime = property(getltctime, "I'm the 'ltctime' property")
302 302 ippSeconds = property(get_ippSeconds, set_ippSeconds)
303 303 dtype = property(get_dtype, set_dtype)
304 304
305 305 class Voltage(JROData):
306 306
307 307 #data es un numpy array de 2 dmensiones (canales, alturas)
308 308 data = None
309 309
310 310 def __init__(self):
311 311 '''
312 312 Constructor
313 313 '''
314 314
315 315 self.radarControllerHeaderObj = RadarControllerHeader()
316 316
317 317 self.systemHeaderObj = SystemHeader()
318 318
319 319 self.type = "Voltage"
320 320
321 321 self.data = None
322 322
323 323 # self.dtype = None
324 324
325 325 # self.nChannels = 0
326 326
327 327 # self.nHeights = 0
328 328
329 329 self.nProfiles = None
330 330
331 331 self.heightList = None
332 332
333 333 self.channelList = None
334 334
335 335 # self.channelIndexList = None
336 336
337 337 self.flagNoData = True
338 338
339 339 self.flagTimeBlock = False
340 340
341 341 self.utctime = None
342 342
343 343 self.timeZone = None
344 344
345 345 self.dstFlag = None
346 346
347 347 self.errorCount = None
348 348
349 349 self.nCohInt = None
350 350
351 351 self.blocksize = None
352 352
353 353 self.flagDecodeData = False #asumo q la data no esta decodificada
354 354
355 355 self.flagDeflipData = False #asumo q la data no esta sin flip
356 356
357 357 self.flagShiftFFT = False
358 358
359 359
360 360 def getNoisebyHildebrand(self):
361 361 """
362 362 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
363 363
364 364 Return:
365 365 noiselevel
366 366 """
367 367
368 368 for channel in range(self.nChannels):
369 369 daux = self.data_spc[channel,:,:]
370 370 self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
371 371
372 372 return self.noise
373 373
374 374 def getNoise(self, type = 1):
375 375
376 376 self.noise = numpy.zeros(self.nChannels)
377 377
378 378 if type == 1:
379 379 noise = self.getNoisebyHildebrand()
380 380
381 381 return 10*numpy.log10(noise)
382 382
383 383 noise = property(getNoise, "I'm the 'nHeights' property.")
384 384
385 385 class Spectra(JROData):
386 386
387 387 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
388 388 data_spc = None
389 389
390 390 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
391 391 data_cspc = None
392 392
393 393 #data es un numpy array de 2 dmensiones (canales, alturas)
394 394 data_dc = None
395 395
396 396 nFFTPoints = None
397 397
398 398 # nPairs = None
399 399
400 400 pairsList = None
401 401
402 402 nIncohInt = None
403 403
404 404 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
405 405
406 406 nCohInt = None #se requiere para determinar el valor de timeInterval
407 407
408 408 ippFactor = None
409 409
410 410 def __init__(self):
411 411 '''
412 412 Constructor
413 413 '''
414 414
415 415 self.radarControllerHeaderObj = RadarControllerHeader()
416 416
417 417 self.systemHeaderObj = SystemHeader()
418 418
419 419 self.type = "Spectra"
420 420
421 421 # self.data = None
422 422
423 423 # self.dtype = None
424 424
425 425 # self.nChannels = 0
426 426
427 427 # self.nHeights = 0
428 428
429 429 self.nProfiles = None
430 430
431 431 self.heightList = None
432 432
433 433 self.channelList = None
434 434
435 435 # self.channelIndexList = None
436 436
437 437 self.pairsList = None
438 438
439 439 self.flagNoData = True
440 440
441 441 self.flagTimeBlock = False
442 442
443 443 self.utctime = None
444 444
445 445 self.nCohInt = None
446 446
447 447 self.nIncohInt = None
448 448
449 449 self.blocksize = None
450 450
451 451 self.nFFTPoints = None
452 452
453 453 self.wavelength = None
454 454
455 455 self.flagDecodeData = False #asumo q la data no esta decodificada
456 456
457 457 self.flagDeflipData = False #asumo q la data no esta sin flip
458 458
459 459 self.flagShiftFFT = False
460 460
461 461 self.ippFactor = 1
462 462
463 463 #self.noise = None
464 464
465 465 self.beacon_heiIndexList = []
466 466
467 467 self.noise_estimation = None
468 468
469 469
470 470 def getNoisebyHildebrand(self):
471 471 """
472 472 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
473 473
474 474 Return:
475 475 noiselevel
476 476 """
477 477
478 478 noise = numpy.zeros(self.nChannels)
479 479 for channel in range(self.nChannels):
480 480 daux = self.data_spc[channel,:,:]
481 481 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
482 482
483 483 return noise
484 484
485 485 def getNoise(self):
486 486 if self.noise_estimation != None:
487 487 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
488 488 else:
489 489 noise = self.getNoisebyHildebrand()
490 490 return noise
491 491
492 492
493 493 def getFreqRange(self, extrapoints=0):
494 494
495 495 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
496 496 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
497 497
498 498 return freqrange
499 499
500 500 def getVelRange(self, extrapoints=0):
501 501
502 502 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
503 503 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
504 504
505 505 return velrange
506 506
507 507 def getNPairs(self):
508 508
509 509 return len(self.pairsList)
510 510
511 511 def getPairsIndexList(self):
512 512
513 513 return range(self.nPairs)
514 514
515 515 def getNormFactor(self):
516 516 pwcode = 1
517 517 if self.flagDecodeData:
518 518 pwcode = numpy.sum(self.code[0]**2)
519 519 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
520 520 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
521 521
522 522 return normFactor
523 523
524 524 def getFlagCspc(self):
525 525
526 526 if self.data_cspc == None:
527 527 return True
528 528
529 529 return False
530 530
531 531 def getFlagDc(self):
532 532
533 533 if self.data_dc == None:
534 534 return True
535 535
536 536 return False
537 537
538 538 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
539 539 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
540 540 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
541 541 flag_cspc = property(getFlagCspc)
542 542 flag_dc = property(getFlagDc)
543 543 noise = property(getNoise, "I'm the 'nHeights' property.")
544 544
545 545 class SpectraHeis(Spectra):
546 546
547 547 data_spc = None
548 548
549 549 data_cspc = None
550 550
551 551 data_dc = None
552 552
553 553 nFFTPoints = None
554 554
555 555 # nPairs = None
556 556
557 557 pairsList = None
558 558
559 559 nIncohInt = None
560 560
561 561 def __init__(self):
562 562
563 563 self.radarControllerHeaderObj = RadarControllerHeader()
564 564
565 565 self.systemHeaderObj = SystemHeader()
566 566
567 567 self.type = "SpectraHeis"
568 568
569 569 # self.dtype = None
570 570
571 571 # self.nChannels = 0
572 572
573 573 # self.nHeights = 0
574 574
575 575 self.nProfiles = None
576 576
577 577 self.heightList = None
578 578
579 579 self.channelList = None
580 580
581 581 # self.channelIndexList = None
582 582
583 583 self.flagNoData = True
584 584
585 585 self.flagTimeBlock = False
586 586
587 587 # self.nPairs = 0
588 588
589 589 self.utctime = None
590 590
591 591 self.blocksize = None
592 592
593 593 def getNormFactor(self):
594 594 pwcode = 1
595 595 if self.flagDecodeData:
596 596 pwcode = numpy.sum(self.code[0]**2)
597 597
598 598 normFactor = self.nIncohInt*self.nCohInt*pwcode
599 599
600 600 return normFactor
601 601
602 602 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
603 603
604 604 class Fits:
605 605
606 606 heightList = None
607 607
608 608 channelList = None
609 609
610 610 flagNoData = True
611 611
612 612 flagTimeBlock = False
613 613
614 614 useLocalTime = False
615 615
616 616 utctime = None
617 617
618 618 timeZone = None
619 619
620 620 # ippSeconds = None
621 621
622 622 timeInterval = None
623 623
624 624 nCohInt = None
625 625
626 626 nIncohInt = None
627 627
628 628 noise = None
629 629
630 630 windowOfFilter = 1
631 631
632 632 #Speed of ligth
633 633 C = 3e8
634 634
635 635 frequency = 49.92e6
636 636
637 637 realtime = False
638 638
639 639
640 640 def __init__(self):
641 641
642 642 self.type = "Fits"
643 643
644 644 self.nProfiles = None
645 645
646 646 self.heightList = None
647 647
648 648 self.channelList = None
649 649
650 650 # self.channelIndexList = None
651 651
652 652 self.flagNoData = True
653 653
654 654 self.utctime = None
655 655
656 656 self.nCohInt = None
657 657
658 658 self.nIncohInt = None
659 659
660 660 self.useLocalTime = True
661 661
662 662 # self.utctime = None
663 663 # self.timeZone = None
664 664 # self.ltctime = None
665 665 # self.timeInterval = None
666 666 # self.header = None
667 667 # self.data_header = None
668 668 # self.data = None
669 669 # self.datatime = None
670 670 # self.flagNoData = False
671 671 # self.expName = ''
672 672 # self.nChannels = None
673 673 # self.nSamples = None
674 674 # self.dataBlocksPerFile = None
675 675 # self.comments = ''
676 676 #
677 677
678 678
679 679 def getltctime(self):
680 680
681 681 if self.useLocalTime:
682 682 return self.utctime - self.timeZone*60
683 683
684 684 return self.utctime
685 685
686 686 def getDatatime(self):
687 687
688 688 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
689 689 return datatime
690 690
691 691 def getTimeRange(self):
692 692
693 693 datatime = []
694 694
695 695 datatime.append(self.ltctime)
696 696 datatime.append(self.ltctime + self.timeInterval)
697 697
698 698 datatime = numpy.array(datatime)
699 699
700 700 return datatime
701 701
702 702 def getHeiRange(self):
703 703
704 704 heis = self.heightList
705 705
706 706 return heis
707 707
708 708 def isEmpty(self):
709 709
710 710 return self.flagNoData
711 711
712 712 def getNHeights(self):
713 713
714 714 return len(self.heightList)
715 715
716 716 def getNChannels(self):
717 717
718 718 return len(self.channelList)
719 719
720 720 def getChannelIndexList(self):
721 721
722 722 return range(self.nChannels)
723 723
724 724 def getNoise(self, type = 1):
725 725
726 726 self.noise = numpy.zeros(self.nChannels)
727 727
728 728 if type == 1:
729 729 noise = self.getNoisebyHildebrand()
730 730
731 731 if type == 2:
732 732 noise = self.getNoisebySort()
733 733
734 734 if type == 3:
735 735 noise = self.getNoisebyWindow()
736 736
737 737 return noise
738 738
739 739 datatime = property(getDatatime, "I'm the 'datatime' property")
740 740 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
741 741 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
742 742 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
743 743 noise = property(getNoise, "I'm the 'nHeights' property.")
744 744 datatime = property(getDatatime, "I'm the 'datatime' property")
745 745 ltctime = property(getltctime, "I'm the 'ltctime' property")
746 746
747 747 class Correlation(JROData):
748 748
749 749 noise = None
750 750
751 751 SNR = None
752 752
753 753 pairsAutoCorr = None #Pairs of Autocorrelation
754 754
755 755 #--------------------------------------------------
756 756
757 757 data_corr = None
758 758
759 759 data_volt = None
760 760
761 761 lagT = None # each element value is a profileIndex
762 762
763 763 lagR = None # each element value is in km
764 764
765 765 pairsList = None
766 766
767 767 calculateVelocity = None
768 768
769 769 nPoints = None
770 770
771 771 nAvg = None
772 772
773 773 bufferSize = None
774 774
775 775 def __init__(self):
776 776 '''
777 777 Constructor
778 778 '''
779 779 self.radarControllerHeaderObj = RadarControllerHeader()
780 780
781 781 self.systemHeaderObj = SystemHeader()
782 782
783 783 self.type = "Correlation"
784 784
785 785 self.data = None
786 786
787 787 self.dtype = None
788 788
789 789 self.nProfiles = None
790 790
791 791 self.heightList = None
792 792
793 793 self.channelList = None
794 794
795 795 self.flagNoData = True
796 796
797 797 self.flagTimeBlock = False
798 798
799 799 self.utctime = None
800 800
801 801 self.timeZone = None
802 802
803 803 self.dstFlag = None
804 804
805 805 self.errorCount = None
806 806
807 807 self.blocksize = None
808 808
809 809 self.flagDecodeData = False #asumo q la data no esta decodificada
810 810
811 811 self.flagDeflipData = False #asumo q la data no esta sin flip
812 812
813 813 self.pairsList = None
814 814
815 815 self.nPoints = None
816 816
817 817 def getLagTRange(self, extrapoints=0):
818 818
819 819 lagTRange = self.lagT
820 820 diff = lagTRange[1] - lagTRange[0]
821 821 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
822 822 lagTRange = numpy.hstack((lagTRange, extra))
823 823
824 824 return lagTRange
825 825
826 826 def getLagRRange(self, extrapoints=0):
827 827
828 828 return self.lagR
829 829
830 830 def getPairsList(self):
831 831
832 832 return self.pairsList
833 833
834 834 def getCalculateVelocity(self):
835 835
836 836 return self.calculateVelocity
837 837
838 838 def getNPoints(self):
839 839
840 840 return self.nPoints
841 841
842 842 def getNAvg(self):
843 843
844 844 return self.nAvg
845 845
846 846 def getBufferSize(self):
847 847
848 848 return self.bufferSize
849 849
850 850 def getPairsAutoCorr(self):
851 851 pairsList = self.pairsList
852 852 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
853 853
854 854 for l in range(len(pairsList)):
855 855 firstChannel = pairsList[l][0]
856 856 secondChannel = pairsList[l][1]
857 857
858 858 #Obteniendo pares de Autocorrelacion
859 859 if firstChannel == secondChannel:
860 860 pairsAutoCorr[firstChannel] = int(l)
861 861
862 862 pairsAutoCorr = pairsAutoCorr.astype(int)
863 863
864 864 return pairsAutoCorr
865 865
866 866 def getNoise(self, mode = 2):
867 867
868 868 indR = numpy.where(self.lagR == 0)[0][0]
869 869 indT = numpy.where(self.lagT == 0)[0][0]
870 870
871 871 jspectra0 = self.data_corr[:,:,indR,:]
872 872 jspectra = copy.copy(jspectra0)
873 873
874 874 num_chan = jspectra.shape[0]
875 875 num_hei = jspectra.shape[2]
876 876
877 877 freq_dc = jspectra.shape[1]/2
878 878 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
879 879
880 880 if ind_vel[0]<0:
881 881 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
882 882
883 883 if mode == 1:
884 884 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
885 885
886 886 if mode == 2:
887 887
888 888 vel = numpy.array([-2,-1,1,2])
889 889 xx = numpy.zeros([4,4])
890 890
891 891 for fil in range(4):
892 892 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
893 893
894 894 xx_inv = numpy.linalg.inv(xx)
895 895 xx_aux = xx_inv[0,:]
896 896
897 897 for ich in range(num_chan):
898 898 yy = jspectra[ich,ind_vel,:]
899 899 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
900 900
901 901 junkid = jspectra[ich,freq_dc,:]<=0
902 902 cjunkid = sum(junkid)
903 903
904 904 if cjunkid.any():
905 905 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
906 906
907 907 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
908 908
909 909 return noise
910 910
911 911 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
912 912 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
913 913 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
914 914 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
915 915 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
916 916
917 917
918 918 class Parameters(JROData):
919 919
920 920 #Information from previous data
921 921
922 922 inputUnit = None #Type of data to be processed
923 923
924 924 operation = None #Type of operation to parametrize
925 925
926 926 normFactor = None #Normalization Factor
927 927
928 928 groupList = None #List of Pairs, Groups, etc
929 929
930 930 #Parameters
931 931
932 932 data_param = None #Parameters obtained
933 933
934 934 data_pre = None #Data Pre Parametrization
935 935
936 data_SNR = None #Signal to Noise Ratio
937
936 938 heightRange = None #Heights
937 939
938 940 abscissaRange = None #Abscissa, can be velocities, lags or time
939 941
940 942 noise = None #Noise Potency
941 943
942 SNR = None #Signal to Noise Ratio
943
944 944 initUtcTime = None #Initial UTC time
945 945
946 946 paramInterval = None #Time interval to calculate Parameters in seconds
947 947
948 948 #Fitting
949 949
950 constants = None
950 data_error = None #Error of the estimation
951 951
952 error = None
952 constants = None
953 953
954 954 library = None
955 955
956 956 #Output signal
957 957
958 958 outputInterval = None #Time interval to calculate output signal in seconds
959 959
960 960 data_output = None #Out signal
961 961
962 962
963 963
964 964 def __init__(self):
965 965 '''
966 966 Constructor
967 967 '''
968 968 self.radarControllerHeaderObj = RadarControllerHeader()
969 969
970 970 self.systemHeaderObj = SystemHeader()
971 971
972 972 self.type = "Parameters"
973 973
974 974 def getTimeRange1(self):
975 975
976 976 datatime = []
977 977
978 978 datatime.append(self.initUtcTime)
979 979 datatime.append(self.initUtcTime + self.outputInterval - 1)
980 980
981 981 datatime = numpy.array(datatime)
982 982
983 983 return datatime
Show file before | Show file after | Hide comments
@@ -1,1178 +1,1178
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from figure import Figure, isRealtime
6 6
7 7 class MomentsPlot(Figure):
8 8
9 9 isConfig = None
10 10 __nsubplots = None
11 11
12 12 WIDTHPROF = None
13 13 HEIGHTPROF = None
14 14 PREFIX = 'prm'
15 15
16 16 def __init__(self):
17 17
18 18 self.isConfig = False
19 19 self.__nsubplots = 1
20 20
21 21 self.WIDTH = 280
22 22 self.HEIGHT = 250
23 23 self.WIDTHPROF = 120
24 24 self.HEIGHTPROF = 0
25 25 self.counter_imagwr = 0
26 26
27 27 self.PLOT_CODE = 1
28 28 self.FTP_WEI = None
29 29 self.EXP_CODE = None
30 30 self.SUB_EXP_CODE = None
31 31 self.PLOT_POS = None
32 32
33 33 def getSubplots(self):
34 34
35 35 ncol = int(numpy.sqrt(self.nplots)+0.9)
36 36 nrow = int(self.nplots*1./ncol + 0.9)
37 37
38 38 return nrow, ncol
39 39
40 40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
41 41
42 42 self.__showprofile = showprofile
43 43 self.nplots = nplots
44 44
45 45 ncolspan = 1
46 46 colspan = 1
47 47 if showprofile:
48 48 ncolspan = 3
49 49 colspan = 2
50 50 self.__nsubplots = 2
51 51
52 52 self.createFigure(id = id,
53 53 wintitle = wintitle,
54 54 widthplot = self.WIDTH + self.WIDTHPROF,
55 55 heightplot = self.HEIGHT + self.HEIGHTPROF,
56 56 show=show)
57 57
58 58 nrow, ncol = self.getSubplots()
59 59
60 60 counter = 0
61 61 for y in range(nrow):
62 62 for x in range(ncol):
63 63
64 64 if counter >= self.nplots:
65 65 break
66 66
67 67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
68 68
69 69 if showprofile:
70 70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
71 71
72 72 counter += 1
73 73
74 74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
75 75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
76 76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
77 77 server=None, folder=None, username=None, password=None,
78 78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
79 79
80 80 """
81 81
82 82 Input:
83 83 dataOut :
84 84 id :
85 85 wintitle :
86 86 channelList :
87 87 showProfile :
88 88 xmin : None,
89 89 xmax : None,
90 90 ymin : None,
91 91 ymax : None,
92 92 zmin : None,
93 93 zmax : None
94 94 """
95 95
96 96 if dataOut.flagNoData:
97 97 return None
98 98
99 99 if realtime:
100 100 if not(isRealtime(utcdatatime = dataOut.utctime)):
101 101 print 'Skipping this plot function'
102 102 return
103 103
104 104 if channelList == None:
105 105 channelIndexList = dataOut.channelIndexList
106 106 else:
107 107 channelIndexList = []
108 108 for channel in channelList:
109 109 if channel not in dataOut.channelList:
110 110 raise ValueError, "Channel %d is not in dataOut.channelList"
111 111 channelIndexList.append(dataOut.channelList.index(channel))
112 112
113 113 factor = dataOut.normFactor
114 114 x = dataOut.abscissaRange
115 115 y = dataOut.heightRange
116 116
117 117 z = dataOut.data_pre[channelIndexList,:,:]/factor
118 118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
119 119 avg = numpy.average(z, axis=1)
120 120 noise = dataOut.noise/factor
121 121
122 122 zdB = 10*numpy.log10(z)
123 123 avgdB = 10*numpy.log10(avg)
124 124 noisedB = 10*numpy.log10(noise)
125 125
126 126 #thisDatetime = dataOut.datatime
127 127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
128 128 title = wintitle + " Parameters"
129 129 xlabel = "Velocity (m/s)"
130 130 ylabel = "Range (Km)"
131 131
132 132 if not self.isConfig:
133 133
134 134 nplots = len(channelIndexList)
135 135
136 136 self.setup(id=id,
137 137 nplots=nplots,
138 138 wintitle=wintitle,
139 139 showprofile=showprofile,
140 140 show=show)
141 141
142 142 if xmin == None: xmin = numpy.nanmin(x)
143 143 if xmax == None: xmax = numpy.nanmax(x)
144 144 if ymin == None: ymin = numpy.nanmin(y)
145 145 if ymax == None: ymax = numpy.nanmax(y)
146 146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
147 147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
148 148
149 149 self.FTP_WEI = ftp_wei
150 150 self.EXP_CODE = exp_code
151 151 self.SUB_EXP_CODE = sub_exp_code
152 152 self.PLOT_POS = plot_pos
153 153
154 154 self.isConfig = True
155 155
156 156 self.setWinTitle(title)
157 157
158 158 for i in range(self.nplots):
159 159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
160 160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
161 161 axes = self.axesList[i*self.__nsubplots]
162 162 axes.pcolor(x, y, zdB[i,:,:],
163 163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
164 164 xlabel=xlabel, ylabel=ylabel, title=title,
165 165 ticksize=9, cblabel='')
166 166 #Mean Line
167 167 mean = dataOut.data_param[i, 1, :]
168 168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
169 169
170 170 if self.__showprofile:
171 171 axes = self.axesList[i*self.__nsubplots +1]
172 172 axes.pline(avgdB[i], y,
173 173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
174 174 xlabel='dB', ylabel='', title='',
175 175 ytick_visible=False,
176 176 grid='x')
177 177
178 178 noiseline = numpy.repeat(noisedB[i], len(y))
179 179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
180 180
181 181 self.draw()
182 182
183 183 if figfile == None:
184 184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
185 185 figfile = self.getFilename(name = str_datetime)
186 186
187 187 if figpath != '':
188 188 self.counter_imagwr += 1
189 189 if (self.counter_imagwr>=wr_period):
190 190 # store png plot to local folder
191 191 self.saveFigure(figpath, figfile)
192 192 # store png plot to FTP server according to RT-Web format
193 193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
194 194 ftp_filename = os.path.join(figpath, name)
195 195 self.saveFigure(figpath, ftp_filename)
196 196 self.counter_imagwr = 0
197 197
198 198 class SkyMapPlot(Figure):
199 199
200 200 __isConfig = None
201 201 __nsubplots = None
202 202
203 203 WIDTHPROF = None
204 204 HEIGHTPROF = None
205 205 PREFIX = 'prm'
206 206
207 207 def __init__(self):
208 208
209 209 self.__isConfig = False
210 210 self.__nsubplots = 1
211 211
212 212 # self.WIDTH = 280
213 213 # self.HEIGHT = 250
214 214 self.WIDTH = 600
215 215 self.HEIGHT = 600
216 216 self.WIDTHPROF = 120
217 217 self.HEIGHTPROF = 0
218 218 self.counter_imagwr = 0
219 219
220 220 self.PLOT_CODE = 1
221 221 self.FTP_WEI = None
222 222 self.EXP_CODE = None
223 223 self.SUB_EXP_CODE = None
224 224 self.PLOT_POS = None
225 225
226 226 def getSubplots(self):
227 227
228 228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 229 nrow = int(self.nplots*1./ncol + 0.9)
230 230
231 231 return nrow, ncol
232 232
233 233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234 234
235 235 self.__showprofile = showprofile
236 236 self.nplots = nplots
237 237
238 238 ncolspan = 1
239 239 colspan = 1
240 240
241 241 self.createFigure(id = id,
242 242 wintitle = wintitle,
243 243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 245 show=show)
246 246
247 247 nrow, ncol = 1,1
248 248 counter = 0
249 249 x = 0
250 250 y = 0
251 251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252 252
253 253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
255 255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 256 server=None, folder=None, username=None, password=None,
257 257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258 258
259 259 """
260 260
261 261 Input:
262 262 dataOut :
263 263 id :
264 264 wintitle :
265 265 channelList :
266 266 showProfile :
267 267 xmin : None,
268 268 xmax : None,
269 269 ymin : None,
270 270 ymax : None,
271 271 zmin : None,
272 272 zmax : None
273 273 """
274 274
275 275 arrayParameters = dataOut.data_param
276 276 error = arrayParameters[:,-1]
277 277 indValid = numpy.where(error == 0)[0]
278 278 finalMeteor = arrayParameters[indValid,:]
279 279 finalAzimuth = finalMeteor[:,4]
280 280 finalZenith = finalMeteor[:,5]
281 281
282 282 x = finalAzimuth*numpy.pi/180
283 283 y = finalZenith
284 284
285 285
286 286 #thisDatetime = dataOut.datatime
287 287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
288 288 title = wintitle + " Parameters"
289 289 xlabel = "Zonal Zenith Angle (deg) "
290 290 ylabel = "Meridional Zenith Angle (deg)"
291 291
292 292 if not self.__isConfig:
293 293
294 294 nplots = 1
295 295
296 296 self.setup(id=id,
297 297 nplots=nplots,
298 298 wintitle=wintitle,
299 299 showprofile=showprofile,
300 300 show=show)
301 301
302 302 self.FTP_WEI = ftp_wei
303 303 self.EXP_CODE = exp_code
304 304 self.SUB_EXP_CODE = sub_exp_code
305 305 self.PLOT_POS = plot_pos
306 306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
307 307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
308 308 self.__isConfig = True
309 309
310 310 self.setWinTitle(title)
311 311
312 312 i = 0
313 313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
314 314
315 315 axes = self.axesList[i*self.__nsubplots]
316 316 nevents = axes.x_buffer.shape[0] + x.shape[0]
317 317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
318 318 axes.polar(x, y,
319 319 title=title, xlabel=xlabel, ylabel=ylabel,
320 320 ticksize=9, cblabel='')
321 321
322 322 self.draw()
323 323
324 324 if save:
325 325
326 326 self.counter_imagwr += 1
327 327 if (self.counter_imagwr==wr_period):
328 328
329 329 if figfile == None:
330 330 figfile = self.getFilename(name = self.name)
331 331 self.saveFigure(figpath, figfile)
332 332
333 333 if ftp:
334 334 #provisionalmente envia archivos en el formato de la web en tiempo real
335 335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
336 336 path = '%s%03d' %(self.PREFIX, self.id)
337 337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
338 338 self.saveFigure(figpath, ftp_file)
339 339 ftp_filename = os.path.join(figpath,ftp_file)
340 340
341 341
342 342 try:
343 343 self.sendByFTP(ftp_filename, server, folder, username, password)
344 344 except:
345 345 self.counter_imagwr = 0
346 346 raise ValueError, 'Error FTP'
347 347
348 348 self.counter_imagwr = 0
349 349
350 350
351 351 class WindProfilerPlot(Figure):
352 352
353 353 __isConfig = None
354 354 __nsubplots = None
355 355
356 356 WIDTHPROF = None
357 357 HEIGHTPROF = None
358 358 PREFIX = 'wind'
359 359
360 360 def __init__(self):
361 361
362 362 self.timerange = 2*60*60
363 363 self.__isConfig = False
364 364 self.__nsubplots = 1
365 365
366 366 self.WIDTH = 800
367 367 self.HEIGHT = 150
368 368 self.WIDTHPROF = 120
369 369 self.HEIGHTPROF = 0
370 370 self.counter_imagwr = 0
371 371
372 372 self.PLOT_CODE = 0
373 373 self.FTP_WEI = None
374 374 self.EXP_CODE = None
375 375 self.SUB_EXP_CODE = None
376 376 self.PLOT_POS = None
377 377 self.tmin = None
378 378 self.tmax = None
379 379
380 380 self.xmin = None
381 381 self.xmax = None
382 382
383 383 self.figfile = None
384 384
385 385 def getSubplots(self):
386 386
387 387 ncol = 1
388 388 nrow = self.nplots
389 389
390 390 return nrow, ncol
391 391
392 392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
393 393
394 394 self.__showprofile = showprofile
395 395 self.nplots = nplots
396 396
397 397 ncolspan = 1
398 398 colspan = 1
399 399
400 400 self.createFigure(id = id,
401 401 wintitle = wintitle,
402 402 widthplot = self.WIDTH + self.WIDTHPROF,
403 403 heightplot = self.HEIGHT + self.HEIGHTPROF,
404 404 show=show)
405 405
406 406 nrow, ncol = self.getSubplots()
407 407
408 408 counter = 0
409 409 for y in range(nrow):
410 410 if counter >= self.nplots:
411 411 break
412 412
413 413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
414 414 counter += 1
415 415
416 416 def run(self, dataOut, id, wintitle="", channelList=None,
417 417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
418 418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
419 419 timerange=None, SNRthresh = None,
420 420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
421 421 server=None, folder=None, username=None, password=None,
422 422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
423 423 """
424 424
425 425 Input:
426 426 dataOut :
427 427 id :
428 428 wintitle :
429 429 channelList :
430 430 showProfile :
431 431 xmin : None,
432 432 xmax : None,
433 433 ymin : None,
434 434 ymax : None,
435 435 zmin : None,
436 436 zmax : None
437 437 """
438 438
439 439 if channelList == None:
440 440 channelIndexList = dataOut.channelIndexList
441 441 else:
442 442 channelIndexList = []
443 443 for channel in channelList:
444 444 if channel not in dataOut.channelList:
445 445 raise ValueError, "Channel %d is not in dataOut.channelList"
446 446 channelIndexList.append(dataOut.channelList.index(channel))
447 447
448 448 if timerange != None:
449 449 self.timerange = timerange
450 450
451 451 tmin = None
452 452 tmax = None
453 453
454 454 x = dataOut.getTimeRange1()
455 455 # y = dataOut.heightRange
456 456 y = dataOut.heightRange
457 457
458 458 z = dataOut.data_output.copy()
459 459 nplots = z.shape[0] #Number of wind dimensions estimated
460 460 nplotsw = nplots
461 461
462 462 #If there is a SNR function defined
463 if dataOut.SNR != None:
463 if dataOut.data_SNR != None:
464 464 nplots += 1
465 SNR = dataOut.SNR
465 SNR = dataOut.data_SNR
466 466 SNRavg = numpy.average(SNR, axis=0)
467 467
468 468 SNRdB = 10*numpy.log10(SNR)
469 469 SNRavgdB = 10*numpy.log10(SNRavg)
470 470
471 471 if SNRthresh == None: SNRthresh = -5.0
472 472 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
473 473
474 474 for i in range(nplotsw):
475 475 z[i,ind] = numpy.nan
476 476
477 477
478 478 showprofile = False
479 479 # thisDatetime = dataOut.datatime
480 480 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
481 481 title = wintitle + "Wind"
482 482 xlabel = ""
483 483 ylabel = "Range (Km)"
484 484
485 485 if not self.__isConfig:
486 486
487 487 self.setup(id=id,
488 488 nplots=nplots,
489 489 wintitle=wintitle,
490 490 showprofile=showprofile,
491 491 show=show)
492 492
493 493 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
494 494
495 495 if ymin == None: ymin = numpy.nanmin(y)
496 496 if ymax == None: ymax = numpy.nanmax(y)
497 497
498 498 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
499 499 #if numpy.isnan(zmax): zmax = 50
500 500 if zmin == None: zmin = -zmax
501 501
502 502 if nplotsw == 3:
503 503 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
504 504 if zmin_ver == None: zmin_ver = -zmax_ver
505 505
506 if dataOut.SNR != None:
506 if dataOut.data_SNR != None:
507 507 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
508 508 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
509 509
510 510 self.FTP_WEI = ftp_wei
511 511 self.EXP_CODE = exp_code
512 512 self.SUB_EXP_CODE = sub_exp_code
513 513 self.PLOT_POS = plot_pos
514 514
515 515 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
516 516 self.__isConfig = True
517 517
518 518
519 519 self.setWinTitle(title)
520 520
521 521 if ((self.xmax - x[1]) < (x[1]-x[0])):
522 522 x[1] = self.xmax
523 523
524 524 strWind = ['Zonal', 'Meridional', 'Vertical']
525 525 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
526 526 zmaxVector = [zmax, zmax, zmax_ver]
527 527 zminVector = [zmin, zmin, zmin_ver]
528 528 windFactor = [1,1,100]
529 529
530 530 for i in range(nplotsw):
531 531
532 532 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
533 533 axes = self.axesList[i*self.__nsubplots]
534 534
535 535 z1 = z[i,:].reshape((1,-1))*windFactor[i]
536 536
537 537 axes.pcolorbuffer(x, y, z1,
538 538 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
539 539 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
540 540 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
541 541
542 if dataOut.SNR != None:
542 if dataOut.data_SNR != None:
543 543 i += 1
544 544 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
545 545 axes = self.axesList[i*self.__nsubplots]
546 546
547 547 SNRavgdB = SNRavgdB.reshape((1,-1))
548 548
549 549 axes.pcolorbuffer(x, y, SNRavgdB,
550 550 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
551 551 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
552 552 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
553 553
554 554 self.draw()
555 555
556 556 if self.figfile == None:
557 557 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
558 558 self.figfile = self.getFilename(name = str_datetime)
559 559
560 560 if figpath != '':
561 561
562 562 self.counter_imagwr += 1
563 563 if (self.counter_imagwr>=wr_period):
564 564 # store png plot to local folder
565 565 self.saveFigure(figpath, self.figfile)
566 566 # store png plot to FTP server according to RT-Web format
567 567 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
568 568 ftp_filename = os.path.join(figpath, name)
569 569 self.saveFigure(figpath, ftp_filename)
570 570
571 571 self.counter_imagwr = 0
572 572
573 573 if x[1] >= self.axesList[0].xmax:
574 574 self.counter_imagwr = wr_period
575 575 self.__isConfig = False
576 576 self.figfile = None
577 577
578 578
579 579 class ParametersPlot(Figure):
580 580
581 581 __isConfig = None
582 582 __nsubplots = None
583 583
584 584 WIDTHPROF = None
585 585 HEIGHTPROF = None
586 586 PREFIX = 'prm'
587 587
588 588 def __init__(self):
589 589
590 590 self.timerange = 2*60*60
591 591 self.__isConfig = False
592 592 self.__nsubplots = 1
593 593
594 594 self.WIDTH = 800
595 595 self.HEIGHT = 150
596 596 self.WIDTHPROF = 120
597 597 self.HEIGHTPROF = 0
598 598 self.counter_imagwr = 0
599 599
600 600 self.PLOT_CODE = 0
601 601 self.FTP_WEI = None
602 602 self.EXP_CODE = None
603 603 self.SUB_EXP_CODE = None
604 604 self.PLOT_POS = None
605 605 self.tmin = None
606 606 self.tmax = None
607 607
608 608 self.xmin = None
609 609 self.xmax = None
610 610
611 611 self.figfile = None
612 612
613 613 def getSubplots(self):
614 614
615 615 ncol = 1
616 616 nrow = self.nplots
617 617
618 618 return nrow, ncol
619 619
620 620 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
621 621
622 622 self.__showprofile = showprofile
623 623 self.nplots = nplots
624 624
625 625 ncolspan = 1
626 626 colspan = 1
627 627
628 628 self.createFigure(id = id,
629 629 wintitle = wintitle,
630 630 widthplot = self.WIDTH + self.WIDTHPROF,
631 631 heightplot = self.HEIGHT + self.HEIGHTPROF,
632 632 show=show)
633 633
634 634 nrow, ncol = self.getSubplots()
635 635
636 636 counter = 0
637 637 for y in range(nrow):
638 638 for x in range(ncol):
639 639
640 640 if counter >= self.nplots:
641 641 break
642 642
643 643 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
644 644
645 645 if showprofile:
646 646 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
647 647
648 648 counter += 1
649 649
650 650 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
651 651 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
652 652 SNRmin = None, SNRmax = None, SNRthresh = None, paramIndex = None, onlyPositive = False,
653 653 zlabel = "", parameterName = "",
654 654 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
655 655 server=None, folder=None, username=None, password=None,
656 656 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
657 657
658 658 """
659 659
660 660 Input:
661 661 dataOut :
662 662 id :
663 663 wintitle :
664 664 channelList :
665 665 showProfile :
666 666 xmin : None,
667 667 xmax : None,
668 668 ymin : None,
669 669 ymax : None,
670 670 zmin : None,
671 671 zmax : None
672 672 """
673 673
674 674 if channelList == None:
675 675 channelIndexList = dataOut.channelIndexList
676 676 else:
677 677 channelIndexList = []
678 678 for channel in channelList:
679 679 if channel not in dataOut.channelList:
680 680 raise ValueError, "Channel %d is not in dataOut.channelList"
681 681 channelIndexList.append(dataOut.channelList.index(channel))
682 682
683 683 if timerange != None:
684 684 self.timerange = timerange
685 685
686 686 #tmin = None
687 687 #tmax = None
688 688 if paramIndex == None:
689 689 paramIndex = 1
690 690 x = dataOut.getTimeRange1()
691 691 y = dataOut.heightRange
692 692 z = dataOut.data_param[channelIndexList,paramIndex,:].copy()
693 693
694 694 zRange = dataOut.abscissaRange
695 695 nplots = z.shape[0] #Number of wind dimensions estimated
696 696 # thisDatetime = dataOut.datatime
697 697 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
698 698 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
699 699 xlabel = ""
700 700 ylabel = "Range (Km)"
701 701
702 702 if onlyPositive:
703 703 colormap = "jet"
704 704 zmin = 0
705 705 else: colormap = "RdBu_r"
706 706
707 707 if not self.__isConfig:
708 708
709 709 self.setup(id=id,
710 710 nplots=nplots,
711 711 wintitle=wintitle,
712 712 showprofile=showprofile,
713 713 show=show)
714 714
715 715 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
716 716
717 717 if ymin == None: ymin = numpy.nanmin(y)
718 718 if ymax == None: ymax = numpy.nanmax(y)
719 719 if zmin == None: zmin = numpy.nanmin(zRange)
720 720 if zmax == None: zmax = numpy.nanmax(zRange)
721 721
722 if dataOut.SNR != None:
722 if dataOut.data_SNR != None:
723 723 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
724 724 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
725 725
726 726 self.FTP_WEI = ftp_wei
727 727 self.EXP_CODE = exp_code
728 728 self.SUB_EXP_CODE = sub_exp_code
729 729 self.PLOT_POS = plot_pos
730 730
731 731 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
732 732 self.__isConfig = True
733 733 self.figfile = figfile
734 734
735 735 self.setWinTitle(title)
736 736
737 737 if ((self.xmax - x[1]) < (x[1]-x[0])):
738 738 x[1] = self.xmax
739 739
740 740 for i in range(nplots):
741 741 title = "%s Channel %d: %s" %(parameterName, dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
742 742
743 743 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
744 744 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
745 745 axes = self.axesList[i*self.__nsubplots]
746 746 z1 = z[i,:].reshape((1,-1))
747 747 axes.pcolorbuffer(x, y, z1,
748 748 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
749 749 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
750 750 ticksize=9, cblabel=zlabel, cbsize="1%")
751 751
752 752 self.draw()
753 753
754 754 if self.figfile == None:
755 755 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
756 756 self.figfile = self.getFilename(name = str_datetime)
757 757
758 758 if figpath != '':
759 759
760 760 self.counter_imagwr += 1
761 761 if (self.counter_imagwr>=wr_period):
762 762 # store png plot to local folder
763 763 self.saveFigure(figpath, self.figfile)
764 764 # store png plot to FTP server according to RT-Web format
765 765 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
766 766 ftp_filename = os.path.join(figpath, name)
767 767 self.saveFigure(figpath, ftp_filename)
768 768
769 769 self.counter_imagwr = 0
770 770
771 771 if x[1] >= self.axesList[0].xmax:
772 772 self.counter_imagwr = wr_period
773 773 self.__isConfig = False
774 774 self.figfile = None
775 775
776 776
777 777 class SpectralFittingPlot(Figure):
778 778
779 779 __isConfig = None
780 780 __nsubplots = None
781 781
782 782 WIDTHPROF = None
783 783 HEIGHTPROF = None
784 784 PREFIX = 'prm'
785 785
786 786
787 787 N = None
788 788 ippSeconds = None
789 789
790 790 def __init__(self):
791 791 self.__isConfig = False
792 792 self.__nsubplots = 1
793 793
794 794 self.WIDTH = 450
795 795 self.HEIGHT = 250
796 796 self.WIDTHPROF = 0
797 797 self.HEIGHTPROF = 0
798 798
799 799 def getSubplots(self):
800 800
801 801 ncol = int(numpy.sqrt(self.nplots)+0.9)
802 802 nrow = int(self.nplots*1./ncol + 0.9)
803 803
804 804 return nrow, ncol
805 805
806 806 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
807 807
808 808 showprofile = False
809 809 self.__showprofile = showprofile
810 810 self.nplots = nplots
811 811
812 812 ncolspan = 5
813 813 colspan = 4
814 814 if showprofile:
815 815 ncolspan = 5
816 816 colspan = 4
817 817 self.__nsubplots = 2
818 818
819 819 self.createFigure(id = id,
820 820 wintitle = wintitle,
821 821 widthplot = self.WIDTH + self.WIDTHPROF,
822 822 heightplot = self.HEIGHT + self.HEIGHTPROF,
823 823 show=show)
824 824
825 825 nrow, ncol = self.getSubplots()
826 826
827 827 counter = 0
828 828 for y in range(nrow):
829 829 for x in range(ncol):
830 830
831 831 if counter >= self.nplots:
832 832 break
833 833
834 834 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
835 835
836 836 if showprofile:
837 837 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
838 838
839 839 counter += 1
840 840
841 841 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
842 842 xmin=None, xmax=None, ymin=None, ymax=None,
843 843 save=False, figpath='./', figfile=None, show=True):
844 844
845 845 """
846 846
847 847 Input:
848 848 dataOut :
849 849 id :
850 850 wintitle :
851 851 channelList :
852 852 showProfile :
853 853 xmin : None,
854 854 xmax : None,
855 855 zmin : None,
856 856 zmax : None
857 857 """
858 858
859 859 if cutHeight==None:
860 860 h=270
861 861 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
862 862 cutHeight = dataOut.heightList[heightindex]
863 863
864 864 factor = dataOut.normFactor
865 865 x = dataOut.abscissaRange[:-1]
866 866 #y = dataOut.getHeiRange()
867 867
868 868 z = dataOut.data_pre[:,:,heightindex]/factor
869 869 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
870 870 avg = numpy.average(z, axis=1)
871 871 listChannels = z.shape[0]
872 872
873 873 #Reconstruct Function
874 874 if fit==True:
875 875 groupArray = dataOut.groupList
876 876 listChannels = groupArray.reshape((groupArray.size))
877 877 listChannels.sort()
878 878 spcFitLine = numpy.zeros(z.shape)
879 879 constants = dataOut.constants
880 880
881 881 nGroups = groupArray.shape[0]
882 882 nChannels = groupArray.shape[1]
883 883 nProfiles = z.shape[1]
884 884
885 885 for f in range(nGroups):
886 886 groupChann = groupArray[f,:]
887 887 p = dataOut.data_param[f,:,heightindex]
888 888 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
889 889 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
890 890 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
891 891 spcFitLine[groupChann,:] = fitLineAux
892 892 # spcFitLine = spcFitLine/factor
893 893
894 894 z = z[listChannels,:]
895 895 spcFitLine = spcFitLine[listChannels,:]
896 896 spcFitLinedB = 10*numpy.log10(spcFitLine)
897 897
898 898 zdB = 10*numpy.log10(z)
899 899 #thisDatetime = dataOut.datatime
900 900 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
901 901 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
902 902 xlabel = "Velocity (m/s)"
903 903 ylabel = "Spectrum"
904 904
905 905 if not self.__isConfig:
906 906
907 907 nplots = listChannels.size
908 908
909 909 self.setup(id=id,
910 910 nplots=nplots,
911 911 wintitle=wintitle,
912 912 showprofile=showprofile,
913 913 show=show)
914 914
915 915 if xmin == None: xmin = numpy.nanmin(x)
916 916 if xmax == None: xmax = numpy.nanmax(x)
917 917 if ymin == None: ymin = numpy.nanmin(zdB)
918 918 if ymax == None: ymax = numpy.nanmax(zdB)+2
919 919
920 920 self.__isConfig = True
921 921
922 922 self.setWinTitle(title)
923 923 for i in range(self.nplots):
924 924 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
925 925 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1)
926 926 axes = self.axesList[i*self.__nsubplots]
927 927 if fit == False:
928 928 axes.pline(x, zdB[i,:],
929 929 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
930 930 xlabel=xlabel, ylabel=ylabel, title=title
931 931 )
932 932 if fit == True:
933 933 fitline=spcFitLinedB[i,:]
934 934 y=numpy.vstack([zdB[i,:],fitline] )
935 935 legendlabels=['Data','Fitting']
936 936 axes.pmultilineyaxis(x, y,
937 937 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
938 938 xlabel=xlabel, ylabel=ylabel, title=title,
939 939 legendlabels=legendlabels, marker=None,
940 940 linestyle='solid', grid='both')
941 941
942 942 self.draw()
943 943
944 944 if save:
945 945 date = thisDatetime.strftime("%Y%m%d_%H%M%S")
946 946 if figfile == None:
947 947 figfile = self.getFilename(name = date)
948 948
949 949 self.saveFigure(figpath, figfile)
950 950
951 951
952 952 class EWDriftsPlot(Figure):
953 953
954 954 __isConfig = None
955 955 __nsubplots = None
956 956
957 957 WIDTHPROF = None
958 958 HEIGHTPROF = None
959 959 PREFIX = 'drift'
960 960
961 961 def __init__(self):
962 962
963 963 self.timerange = 2*60*60
964 964 self.isConfig = False
965 965 self.__nsubplots = 1
966 966
967 967 self.WIDTH = 800
968 968 self.HEIGHT = 150
969 969 self.WIDTHPROF = 120
970 970 self.HEIGHTPROF = 0
971 971 self.counter_imagwr = 0
972 972
973 973 self.PLOT_CODE = 0
974 974 self.FTP_WEI = None
975 975 self.EXP_CODE = None
976 976 self.SUB_EXP_CODE = None
977 977 self.PLOT_POS = None
978 978 self.tmin = None
979 979 self.tmax = None
980 980
981 981 self.xmin = None
982 982 self.xmax = None
983 983
984 984 self.figfile = None
985 985
986 986 def getSubplots(self):
987 987
988 988 ncol = 1
989 989 nrow = self.nplots
990 990
991 991 return nrow, ncol
992 992
993 993 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
994 994
995 995 self.__showprofile = showprofile
996 996 self.nplots = nplots
997 997
998 998 ncolspan = 1
999 999 colspan = 1
1000 1000
1001 1001 self.createFigure(id = id,
1002 1002 wintitle = wintitle,
1003 1003 widthplot = self.WIDTH + self.WIDTHPROF,
1004 1004 heightplot = self.HEIGHT + self.HEIGHTPROF,
1005 1005 show=show)
1006 1006
1007 1007 nrow, ncol = self.getSubplots()
1008 1008
1009 1009 counter = 0
1010 1010 for y in range(nrow):
1011 1011 if counter >= self.nplots:
1012 1012 break
1013 1013
1014 1014 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1015 1015 counter += 1
1016 1016
1017 1017 def run(self, dataOut, id, wintitle="", channelList=None,
1018 1018 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1019 1019 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1020 1020 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1021 1021 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1022 1022 server=None, folder=None, username=None, password=None,
1023 1023 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1024 1024 """
1025 1025
1026 1026 Input:
1027 1027 dataOut :
1028 1028 id :
1029 1029 wintitle :
1030 1030 channelList :
1031 1031 showProfile :
1032 1032 xmin : None,
1033 1033 xmax : None,
1034 1034 ymin : None,
1035 1035 ymax : None,
1036 1036 zmin : None,
1037 1037 zmax : None
1038 1038 """
1039 1039
1040 1040 if channelList == None:
1041 1041 channelIndexList = dataOut.channelIndexList
1042 1042 else:
1043 1043 channelIndexList = []
1044 1044 for channel in channelList:
1045 1045 if channel not in dataOut.channelList:
1046 1046 raise ValueError, "Channel %d is not in dataOut.channelList"
1047 1047 channelIndexList.append(dataOut.channelList.index(channel))
1048 1048
1049 1049 if timerange != None:
1050 1050 self.timerange = timerange
1051 1051
1052 1052 tmin = None
1053 1053 tmax = None
1054 1054
1055 1055 x = dataOut.getTimeRange1()
1056 1056 # y = dataOut.heightRange
1057 1057 y = dataOut.heightList
1058 1058
1059 1059 z = dataOut.data_output
1060 1060 nplots = z.shape[0] #Number of wind dimensions estimated
1061 1061 nplotsw = nplots
1062 1062
1063 1063 #If there is a SNR function defined
1064 if dataOut.SNR != None:
1064 if dataOut.data_SNR != None:
1065 1065 nplots += 1
1066 SNR = dataOut.SNR
1066 SNR = dataOut.data_SNR
1067 1067
1068 1068 if SNR_1:
1069 1069 SNR += 1
1070 1070
1071 1071 SNRavg = numpy.average(SNR, axis=0)
1072 1072
1073 1073 SNRdB = 10*numpy.log10(SNR)
1074 1074 SNRavgdB = 10*numpy.log10(SNRavg)
1075 1075
1076 1076 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1077 1077
1078 1078 for i in range(nplotsw):
1079 1079 z[i,ind] = numpy.nan
1080 1080
1081 1081
1082 1082 showprofile = False
1083 1083 # thisDatetime = dataOut.datatime
1084 1084 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1085 1085 title = wintitle + " EW Drifts"
1086 1086 xlabel = ""
1087 1087 ylabel = "Height (Km)"
1088 1088
1089 1089 if not self.__isConfig:
1090 1090
1091 1091 self.setup(id=id,
1092 1092 nplots=nplots,
1093 1093 wintitle=wintitle,
1094 1094 showprofile=showprofile,
1095 1095 show=show)
1096 1096
1097 1097 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1098 1098
1099 1099 if ymin == None: ymin = numpy.nanmin(y)
1100 1100 if ymax == None: ymax = numpy.nanmax(y)
1101 1101
1102 1102 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1103 1103 if zminZonal == None: zminZonal = -zmaxZonal
1104 1104 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1105 1105 if zminVertical == None: zminVertical = -zmaxVertical
1106 1106
1107 if dataOut.SNR != None:
1107 if dataOut.data_SNR != None:
1108 1108 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1109 1109 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1110 1110
1111 1111 self.FTP_WEI = ftp_wei
1112 1112 self.EXP_CODE = exp_code
1113 1113 self.SUB_EXP_CODE = sub_exp_code
1114 1114 self.PLOT_POS = plot_pos
1115 1115
1116 1116 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1117 1117 self.__isConfig = True
1118 1118
1119 1119
1120 1120 self.setWinTitle(title)
1121 1121
1122 1122 if ((self.xmax - x[1]) < (x[1]-x[0])):
1123 1123 x[1] = self.xmax
1124 1124
1125 1125 strWind = ['Zonal','Vertical']
1126 1126 strCb = 'Velocity (m/s)'
1127 1127 zmaxVector = [zmaxZonal, zmaxVertical]
1128 1128 zminVector = [zminZonal, zminVertical]
1129 1129
1130 1130 for i in range(nplotsw):
1131 1131
1132 1132 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1133 1133 axes = self.axesList[i*self.__nsubplots]
1134 1134
1135 1135 z1 = z[i,:].reshape((1,-1))
1136 1136
1137 1137 axes.pcolorbuffer(x, y, z1,
1138 1138 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1139 1139 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1140 1140 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1141 1141
1142 if dataOut.SNR != None:
1142 if dataOut.data_SNR != None:
1143 1143 i += 1
1144 1144 if SNR_1:
1145 1145 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1146 1146 else:
1147 1147 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1148 1148 axes = self.axesList[i*self.__nsubplots]
1149 1149 SNRavgdB = SNRavgdB.reshape((1,-1))
1150 1150
1151 1151 axes.pcolorbuffer(x, y, SNRavgdB,
1152 1152 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1153 1153 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1154 1154 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1155 1155
1156 1156 self.draw()
1157 1157
1158 1158 if self.figfile == None:
1159 1159 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1160 1160 self.figfile = self.getFilename(name = str_datetime)
1161 1161
1162 1162 if figpath != '':
1163 1163
1164 1164 self.counter_imagwr += 1
1165 1165 if (self.counter_imagwr>=wr_period):
1166 1166 # store png plot to local folder
1167 1167 self.saveFigure(figpath, self.figfile)
1168 1168 # store png plot to FTP server according to RT-Web format
1169 1169 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1170 1170 ftp_filename = os.path.join(figpath, name)
1171 1171 self.saveFigure(figpath, ftp_filename)
1172 1172
1173 1173 self.counter_imagwr = 0
1174 1174
1175 1175 if x[1] >= self.axesList[0].xmax:
1176 1176 self.counter_imagwr = wr_period
1177 1177 self.__isConfig = False
1178 1178 self.figfile = None No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,4 +1,5
1 1 from jroIO_voltage import *
2 2 from jroIO_spectra import *
3 3 from jroIO_heispectra import *
4 4 from jroIO_amisr import *
5 from jroIO_HDF5 import * No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1749 +1,1754
1 1 import numpy
2 2 import math
3 3 from scipy import optimize
4 4 from scipy import interpolate
5 5 from scipy import signal
6 6 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10 10 import sys
11 11 import importlib
12 12 import itertools
13 13
14 14 from jroproc_base import ProcessingUnit, Operation
15 15 from model.data.jrodata import Parameters
16 16
17 17
18 18 class ParametersProc(ProcessingUnit):
19 19
20 20 nSeconds = None
21 21
22 22 def __init__(self):
23 23 ProcessingUnit.__init__(self)
24 24
25 25 self.objectDict = {}
26 26 self.buffer = None
27 27 self.firstdatatime = None
28 28 self.profIndex = 0
29 29 self.dataOut = Parameters()
30 30
31 31 def __updateObjFromInput(self):
32 32
33 33 self.dataOut.inputUnit = self.dataIn.type
34 34
35 35 self.dataOut.timeZone = self.dataIn.timeZone
36 36 self.dataOut.dstFlag = self.dataIn.dstFlag
37 37 self.dataOut.errorCount = self.dataIn.errorCount
38 38 self.dataOut.useLocalTime = self.dataIn.useLocalTime
39 39
40 40 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
41 41 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
42 42 self.dataOut.channelList = self.dataIn.channelList
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
45 45 # self.dataOut.nHeights = self.dataIn.nHeights
46 46 # self.dataOut.nChannels = self.dataIn.nChannels
47 47 self.dataOut.nBaud = self.dataIn.nBaud
48 48 self.dataOut.nCode = self.dataIn.nCode
49 49 self.dataOut.code = self.dataIn.code
50 50 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 51 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
52 52 self.dataOut.utctime = self.firstdatatime
53 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
54 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
55 55 # self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 # self.dataOut.nIncohInt = 1
57 57 self.dataOut.ippSeconds = self.dataIn.ippSeconds
58 58 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59 self.dataOut.timeInterval = self.dataIn.timeInterval
60 60 self.dataOut.heightRange = self.dataIn.getHeiRange()
61 61 self.dataOut.frequency = self.dataIn.frequency
62 62
63 63 def run(self, nSeconds = None, nProfiles = None):
64 64
65 65
66 66
67 67 if self.firstdatatime == None:
68 68 self.firstdatatime = self.dataIn.utctime
69 69
70 70 #---------------------- Voltage Data ---------------------------
71 71
72 72 if self.dataIn.type == "Voltage":
73 73 self.dataOut.flagNoData = True
74 74 if nSeconds != None:
75 75 self.nSeconds = nSeconds
76 76 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
77 77
78 78 if self.buffer == None:
79 79 self.buffer = numpy.zeros((self.dataIn.nChannels,
80 80 self.nProfiles,
81 81 self.dataIn.nHeights),
82 82 dtype='complex')
83 83
84 84 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
85 85 self.profIndex += 1
86 86
87 87 if self.profIndex == self.nProfiles:
88 88
89 89 self.__updateObjFromInput()
90 90 self.dataOut.data_pre = self.buffer.copy()
91 91 self.dataOut.paramInterval = nSeconds
92 92 self.dataOut.flagNoData = False
93 93
94 94 self.buffer = None
95 95 self.firstdatatime = None
96 96 self.profIndex = 0
97 97 return
98 98
99 99 #---------------------- Spectra Data ---------------------------
100 100
101 101 if self.dataIn.type == "Spectra":
102 102 self.dataOut.data_pre = self.dataIn.data_spc.copy()
103 103 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
104 104 self.dataOut.noise = self.dataIn.getNoise()
105 105 self.dataOut.normFactor = self.dataIn.normFactor
106 106 self.dataOut.flagNoData = False
107 107
108 108 #---------------------- Correlation Data ---------------------------
109 109
110 110 if self.dataIn.type == "Correlation":
111 111 lagRRange = self.dataIn.lagR
112 112 indR = numpy.where(lagRRange == 0)[0][0]
113 113
114 114 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
115 115 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
116 116 self.dataOut.noise = self.dataIn.noise
117 117 self.dataOut.normFactor = self.dataIn.normFactor
118 self.dataOut.SNR = self.dataIn.SNR
118 self.dataOut.data_SNR = self.dataIn.SNR
119 119 self.dataOut.groupList = self.dataIn.pairsList
120 120 self.dataOut.flagNoData = False
121 121
122 122
123 123 self.__updateObjFromInput()
124 124 self.firstdatatime = None
125 125 self.dataOut.initUtcTime = self.dataIn.ltctime
126 126 self.dataOut.outputInterval = self.dataIn.timeInterval
127 127
128 128 #------------------- Get Moments ----------------------------------
129 129 def GetMoments(self, channelList = None):
130 130 '''
131 131 Function GetMoments()
132 132
133 133 Input:
134 134 channelList : simple channel list to select e.g. [2,3,7]
135 135 self.dataOut.data_pre
136 136 self.dataOut.abscissaRange
137 137 self.dataOut.noise
138 138
139 139 Affected:
140 140 self.dataOut.data_param
141 self.dataOut.SNR
141 self.dataOut.data_SNR
142 142
143 143 '''
144 144 data = self.dataOut.data_pre
145 145 absc = self.dataOut.abscissaRange[:-1]
146 146 noise = self.dataOut.noise
147 147
148 148 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
149 149
150 150 if channelList== None:
151 151 channelList = self.dataIn.channelList
152 152 self.dataOut.channelList = channelList
153 153
154 154 for ind in channelList:
155 155 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
156 156
157 157 self.dataOut.data_param = data_param[:,1:,:]
158 self.dataOut.SNR = data_param[:,0]
158 self.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 163 if (nicoh == None): nicoh = 1
164 164 if (graph == None): graph = 0
165 165 if (smooth == None): smooth = 0
166 166 elif (self.smooth < 3): smooth = 0
167 167
168 168 if (type1 == None): type1 = 0
169 169 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
170 170 if (snrth == None): snrth = -3
171 171 if (dc == None): dc = 0
172 172 if (aliasing == None): aliasing = 0
173 173 if (oldfd == None): oldfd = 0
174 174 if (wwauto == 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 Lags ----------------------------------
234 234
235 235 def GetLags(self):
236 236 '''
237 237 Function GetMoments()
238 238
239 239 Input:
240 240 self.dataOut.data_pre
241 241 self.dataOut.abscissaRange
242 242 self.dataOut.noise
243 243 self.dataOut.normFactor
244 self.dataOut.SNR
244 self.dataOut.data_SNR
245 245 self.dataOut.groupList
246 246 self.dataOut.nChannels
247 247
248 248 Affected:
249 249 self.dataOut.data_param
250 250
251 251 '''
252 252 data = self.dataOut.data_pre
253 253 normFactor = self.dataOut.normFactor
254 254 nHeights = self.dataOut.nHeights
255 255 absc = self.dataOut.abscissaRange[:-1]
256 256 noise = self.dataOut.noise
257 SNR = self.dataOut.SNR
257 SNR = self.dataOut.data_SNR
258 258 pairsList = self.dataOut.groupList
259 259 nChannels = self.dataOut.nChannels
260 260 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
261 261 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
262 262
263 263 dataNorm = numpy.abs(data)
264 264 for l in range(len(pairsList)):
265 265 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
266 266
267 267 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
268 268 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
269 269 return
270 270
271 271 def __getPairsAutoCorr(self, pairsList, nChannels):
272 272
273 273 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
274 274
275 275 for l in range(len(pairsList)):
276 276 firstChannel = pairsList[l][0]
277 277 secondChannel = pairsList[l][1]
278 278
279 279 #Obteniendo pares de Autocorrelacion
280 280 if firstChannel == secondChannel:
281 281 pairsAutoCorr[firstChannel] = int(l)
282 282
283 283 pairsAutoCorr = pairsAutoCorr.astype(int)
284 284
285 285 pairsCrossCorr = range(len(pairsList))
286 286 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
287 287
288 288 return pairsAutoCorr, pairsCrossCorr
289 289
290 290 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
291 291
292 292 Pt0 = data.shape[1]/2
293 293 #Funcion de Autocorrelacion
294 294 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
295 295
296 296 #Obtencion Indice de TauCross
297 297 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
298 298 #Obtencion Indice de TauAuto
299 299 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
300 300 CCValue = data[pairsCrossCorr,Pt0,:]
301 301 for i in range(pairsCrossCorr.size):
302 302 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
303 303
304 304 #Obtencion de TauCross y TauAuto
305 305 tauCross = lagTRange[indCross]
306 306 tauAuto = lagTRange[indAuto]
307 307
308 308 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
309 309
310 310 tauCross[Nan1,Nan2] = numpy.nan
311 311 tauAuto[Nan1,Nan2] = numpy.nan
312 312 tau = numpy.vstack((tauCross,tauAuto))
313 313
314 314 return tau
315 315
316 316 def __calculateLag1Phase(self, data, pairs, lagTRange):
317 317 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
318 318 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
319 319
320 320 phase = numpy.angle(data1[lag1,:])
321 321
322 322 return phase
323 323 #------------------- Detect Meteors ------------------------------
324 324
325 325 def DetectMeteors(self, hei_ref = None, tauindex = 0,
326 326 predefinedPhaseShifts = None, centerReceiverIndex = 2,
327 327 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
328 328 noise_timeStep = 4, noise_multiple = 4,
329 329 multDet_timeLimit = 1, multDet_rangeLimit = 3,
330 330 phaseThresh = 20, SNRThresh = 8,
331 331 hmin = 70, hmax=110, azimuth = 0) :
332 332
333 333 '''
334 334 Function DetectMeteors()
335 335 Project developed with paper:
336 336 HOLDSWORTH ET AL. 2004
337 337
338 338 Input:
339 339 self.dataOut.data_pre
340 340
341 341 centerReceiverIndex: From the channels, which is the center receiver
342 342
343 343 hei_ref: Height reference for the Beacon signal extraction
344 344 tauindex:
345 345 predefinedPhaseShifts: Predefined phase offset for the voltge signals
346 346
347 347 cohDetection: Whether to user Coherent detection or not
348 348 cohDet_timeStep: Coherent Detection calculation time step
349 349 cohDet_thresh: Coherent Detection phase threshold to correct phases
350 350
351 351 noise_timeStep: Noise calculation time step
352 352 noise_multiple: Noise multiple to define signal threshold
353 353
354 354 multDet_timeLimit: Multiple Detection Removal time limit in seconds
355 355 multDet_rangeLimit: Multiple Detection Removal range limit in km
356 356
357 357 phaseThresh: Maximum phase difference between receiver to be consider a meteor
358 358 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
359 359
360 360 hmin: Minimum Height of the meteor to use it in the further wind estimations
361 361 hmax: Maximum Height of the meteor to use it in the further wind estimations
362 362 azimuth: Azimuth angle correction
363 363
364 364 Affected:
365 365 self.dataOut.data_param
366 366
367 367 Rejection Criteria (Errors):
368 368 0: No error; analysis OK
369 369 1: SNR < SNR threshold
370 370 2: angle of arrival (AOA) ambiguously determined
371 371 3: AOA estimate not feasible
372 372 4: Large difference in AOAs obtained from different antenna baselines
373 373 5: echo at start or end of time series
374 374 6: echo less than 5 examples long; too short for analysis
375 375 7: echo rise exceeds 0.3s
376 376 8: echo decay time less than twice rise time
377 377 9: large power level before echo
378 378 10: large power level after echo
379 379 11: poor fit to amplitude for estimation of decay time
380 380 12: poor fit to CCF phase variation for estimation of radial drift velocity
381 381 13: height unresolvable echo: not valid height within 70 to 110 km
382 382 14: height ambiguous echo: more then one possible height within 70 to 110 km
383 383 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
384 384 16: oscilatory echo, indicating event most likely not an underdense echo
385 385
386 386 17: phase difference in meteor Reestimation
387 387
388 388 Data Storage:
389 389 Meteors for Wind Estimation (8):
390 390 Day Hour | Range Height
391 391 Azimuth Zenith errorCosDir
392 392 VelRad errorVelRad
393 393 TypeError
394 394
395 395 '''
396 396 #Get Beacon signal
397 397 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
398 398
399 399 if hei_ref != None:
400 400 newheis = numpy.where(self.dataOut.heightList>hei_ref)
401 401
402 402 heiRang = self.dataOut.getHeiRange()
403 403 #Pairs List
404 404 pairslist = []
405 405 nChannel = self.dataOut.nChannels
406 406 for i in range(nChannel):
407 407 if i != centerReceiverIndex:
408 408 pairslist.append((centerReceiverIndex,i))
409 409
410 410 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
411 411 # see if the user put in pre defined phase shifts
412 412 voltsPShift = self.dataOut.data_pre.copy()
413 413
414 414 if predefinedPhaseShifts != None:
415 415 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
416 416 else:
417 417 #get hardware phase shifts using beacon signal
418 418 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
419 419 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
420 420
421 421 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
422 422 for i in range(self.dataOut.data_pre.shape[0]):
423 423 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
424 424 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
425 425
426 426 #Remove DC
427 427 voltsDC = numpy.mean(voltsPShift,1)
428 428 voltsDC = numpy.mean(voltsDC,1)
429 429 for i in range(voltsDC.shape[0]):
430 430 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
431 431
432 432 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
433 433 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
434 434
435 435 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
436 436 #Coherent Detection
437 437 if cohDetection:
438 438 #use coherent detection to get the net power
439 439 cohDet_thresh = cohDet_thresh*numpy.pi/180
440 440 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
441 441
442 442 #Non-coherent detection!
443 443 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
444 444 #********** END OF COH/NON-COH POWER CALCULATION**********************
445 445
446 446 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
447 447 #Get noise
448 448 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
449 449 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
450 450 #Get signal threshold
451 451 signalThresh = noise_multiple*noise
452 452 #Meteor echoes detection
453 453 listMeteors = self.__findMeteors(powerNet, signalThresh)
454 454 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
455 455
456 456 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
457 457 #Parameters
458 458 heiRange = self.dataOut.getHeiRange()
459 459 rangeInterval = heiRange[1] - heiRange[0]
460 460 rangeLimit = multDet_rangeLimit/rangeInterval
461 461 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
462 462 #Multiple detection removals
463 463 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
464 464 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
465 465
466 466 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
467 467 #Parameters
468 468 phaseThresh = phaseThresh*numpy.pi/180
469 469 thresh = [phaseThresh, noise_multiple, SNRThresh]
470 470 #Meteor reestimation (Errors N 1, 6, 12, 17)
471 471 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
472 472 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
473 473 #Estimation of decay times (Errors N 7, 8, 11)
474 474 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
475 475 #******************* END OF METEOR REESTIMATION *******************
476 476
477 477 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
478 478 #Calculating Radial Velocity (Error N 15)
479 479 radialStdThresh = 10
480 480 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
481 481
482 482 if len(listMeteors4) > 0:
483 483 #Setting New Array
484 484 date = repr(self.dataOut.datatime)
485 485 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
486 486
487 487 #Calculate AOA (Error N 3, 4)
488 488 #JONES ET AL. 1998
489 489 AOAthresh = numpy.pi/8
490 490 error = arrayParameters[:,-1]
491 491 phases = -arrayMeteors4[:,9:13]
492 492 pairsList = []
493 493 pairsList.append((0,3))
494 494 pairsList.append((1,2))
495 495 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
496 496
497 497 #Calculate Heights (Error N 13 and 14)
498 498 error = arrayParameters[:,-1]
499 499 Ranges = arrayParameters[:,2]
500 500 zenith = arrayParameters[:,5]
501 501 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
502 502 #********************* END OF PARAMETERS CALCULATION **************************
503 503
504 504 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
505 505 self.dataOut.data_param = arrayParameters
506 506
507 507 return
508 508
509 509 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
510 510
511 511 minIndex = min(newheis[0])
512 512 maxIndex = max(newheis[0])
513 513
514 514 voltage = voltage0[:,:,minIndex:maxIndex+1]
515 515 nLength = voltage.shape[1]/n
516 516 nMin = 0
517 517 nMax = 0
518 518 phaseOffset = numpy.zeros((len(pairslist),n))
519 519
520 520 for i in range(n):
521 521 nMax += nLength
522 522 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
523 523 phaseCCF = numpy.mean(phaseCCF, axis = 2)
524 524 phaseOffset[:,i] = phaseCCF.transpose()
525 525 nMin = nMax
526 526 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
527 527
528 528 #Remove Outliers
529 529 factor = 2
530 530 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
531 531 dw = numpy.std(wt,axis = 1)
532 532 dw = dw.reshape((dw.size,1))
533 533 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
534 534 phaseOffset[ind] = numpy.nan
535 535 phaseOffset = stats.nanmean(phaseOffset, axis=1)
536 536
537 537 return phaseOffset
538 538
539 539 def __shiftPhase(self, data, phaseShift):
540 540 #this will shift the phase of a complex number
541 541 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
542 542 return dataShifted
543 543
544 544 def __estimatePhaseDifference(self, array, pairslist):
545 545 nChannel = array.shape[0]
546 546 nHeights = array.shape[2]
547 547 numPairs = len(pairslist)
548 548 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
549 549 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
550 550
551 551 #Correct phases
552 552 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
553 553 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
554 554
555 555 if indDer[0].shape[0] > 0:
556 556 for i in range(indDer[0].shape[0]):
557 557 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
558 558 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
559 559
560 560 # for j in range(numSides):
561 561 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
562 562 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
563 563 #
564 564 #Linear
565 565 phaseInt = numpy.zeros((numPairs,1))
566 566 angAllCCF = phaseCCF[:,[0,1,3,4],0]
567 567 for j in range(numPairs):
568 568 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
569 569 phaseInt[j] = fit[1]
570 570 #Phase Differences
571 571 phaseDiff = phaseInt - phaseCCF[:,2,:]
572 572 phaseArrival = phaseInt.reshape(phaseInt.size)
573 573
574 574 #Dealias
575 575 indAlias = numpy.where(phaseArrival > numpy.pi)
576 576 phaseArrival[indAlias] -= 2*numpy.pi
577 577 indAlias = numpy.where(phaseArrival < -numpy.pi)
578 578 phaseArrival[indAlias] += 2*numpy.pi
579 579
580 580 return phaseDiff, phaseArrival
581 581
582 582 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
583 583 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
584 584 #find the phase shifts of each channel over 1 second intervals
585 585 #only look at ranges below the beacon signal
586 586 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
587 587 numBlocks = int(volts.shape[1]/numProfPerBlock)
588 588 numHeights = volts.shape[2]
589 589 nChannel = volts.shape[0]
590 590 voltsCohDet = volts.copy()
591 591
592 592 pairsarray = numpy.array(pairslist)
593 593 indSides = pairsarray[:,1]
594 594 # indSides = numpy.array(range(nChannel))
595 595 # indSides = numpy.delete(indSides, indCenter)
596 596 #
597 597 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
598 598 listBlocks = numpy.array_split(volts, numBlocks, 1)
599 599
600 600 startInd = 0
601 601 endInd = 0
602 602
603 603 for i in range(numBlocks):
604 604 startInd = endInd
605 605 endInd = endInd + listBlocks[i].shape[1]
606 606
607 607 arrayBlock = listBlocks[i]
608 608 # arrayBlockCenter = listCenter[i]
609 609
610 610 #Estimate the Phase Difference
611 611 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
612 612 #Phase Difference RMS
613 613 arrayPhaseRMS = numpy.abs(phaseDiff)
614 614 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
615 615 indPhase = numpy.where(phaseRMSaux==4)
616 616 #Shifting
617 617 if indPhase[0].shape[0] > 0:
618 618 for j in range(indSides.size):
619 619 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
620 620 voltsCohDet[:,startInd:endInd,:] = arrayBlock
621 621
622 622 return voltsCohDet
623 623
624 624 def __calculateCCF(self, volts, pairslist ,laglist):
625 625
626 626 nHeights = volts.shape[2]
627 627 nPoints = volts.shape[1]
628 628 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
629 629
630 630 for i in range(len(pairslist)):
631 631 volts1 = volts[pairslist[i][0]]
632 632 volts2 = volts[pairslist[i][1]]
633 633
634 634 for t in range(len(laglist)):
635 635 idxT = laglist[t]
636 636 if idxT >= 0:
637 637 vStacked = numpy.vstack((volts2[idxT:,:],
638 638 numpy.zeros((idxT, nHeights),dtype='complex')))
639 639 else:
640 640 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
641 641 volts2[:(nPoints + idxT),:]))
642 642 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
643 643
644 644 vStacked = None
645 645 return voltsCCF
646 646
647 647 def __getNoise(self, power, timeSegment, timeInterval):
648 648 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
649 649 numBlocks = int(power.shape[0]/numProfPerBlock)
650 650 numHeights = power.shape[1]
651 651
652 652 listPower = numpy.array_split(power, numBlocks, 0)
653 653 noise = numpy.zeros((power.shape[0], power.shape[1]))
654 654 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
655 655
656 656 startInd = 0
657 657 endInd = 0
658 658
659 659 for i in range(numBlocks): #split por canal
660 660 startInd = endInd
661 661 endInd = endInd + listPower[i].shape[0]
662 662
663 663 arrayBlock = listPower[i]
664 664 noiseAux = numpy.mean(arrayBlock, 0)
665 665 # noiseAux = numpy.median(noiseAux)
666 666 # noiseAux = numpy.mean(arrayBlock)
667 667 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
668 668
669 669 noiseAux1 = numpy.mean(arrayBlock)
670 670 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
671 671
672 672 return noise, noise1
673 673
674 674 def __findMeteors(self, power, thresh):
675 675 nProf = power.shape[0]
676 676 nHeights = power.shape[1]
677 677 listMeteors = []
678 678
679 679 for i in range(nHeights):
680 680 powerAux = power[:,i]
681 681 threshAux = thresh[:,i]
682 682
683 683 indUPthresh = numpy.where(powerAux > threshAux)[0]
684 684 indDNthresh = numpy.where(powerAux <= threshAux)[0]
685 685
686 686 j = 0
687 687
688 688 while (j < indUPthresh.size - 2):
689 689 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
690 690 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
691 691 indDNthresh = indDNthresh[indDNAux]
692 692
693 693 if (indDNthresh.size > 0):
694 694 indEnd = indDNthresh[0] - 1
695 695 indInit = indUPthresh[j]
696 696
697 697 meteor = powerAux[indInit:indEnd + 1]
698 698 indPeak = meteor.argmax() + indInit
699 699 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
700 700
701 701 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
702 702 j = numpy.where(indUPthresh == indEnd)[0] + 1
703 703 else: j+=1
704 704 else: j+=1
705 705
706 706 return listMeteors
707 707
708 708 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
709 709
710 710 arrayMeteors = numpy.asarray(listMeteors)
711 711 listMeteors1 = []
712 712
713 713 while arrayMeteors.shape[0] > 0:
714 714 FLAs = arrayMeteors[:,4]
715 715 maxFLA = FLAs.argmax()
716 716 listMeteors1.append(arrayMeteors[maxFLA,:])
717 717
718 718 MeteorInitTime = arrayMeteors[maxFLA,1]
719 719 MeteorEndTime = arrayMeteors[maxFLA,3]
720 720 MeteorHeight = arrayMeteors[maxFLA,0]
721 721
722 722 #Check neighborhood
723 723 maxHeightIndex = MeteorHeight + rangeLimit
724 724 minHeightIndex = MeteorHeight - rangeLimit
725 725 minTimeIndex = MeteorInitTime - timeLimit
726 726 maxTimeIndex = MeteorEndTime + timeLimit
727 727
728 728 #Check Heights
729 729 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
730 730 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
731 731 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
732 732
733 733 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
734 734
735 735 return listMeteors1
736 736
737 737 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
738 738 numHeights = volts.shape[2]
739 739 nChannel = volts.shape[0]
740 740
741 741 thresholdPhase = thresh[0]
742 742 thresholdNoise = thresh[1]
743 743 thresholdDB = float(thresh[2])
744 744
745 745 thresholdDB1 = 10**(thresholdDB/10)
746 746 pairsarray = numpy.array(pairslist)
747 747 indSides = pairsarray[:,1]
748 748
749 749 pairslist1 = list(pairslist)
750 750 pairslist1.append((0,1))
751 751 pairslist1.append((3,4))
752 752
753 753 listMeteors1 = []
754 754 listPowerSeries = []
755 755 listVoltageSeries = []
756 756 #volts has the war data
757 757
758 758 if frequency == 30e6:
759 759 timeLag = 45*10**-3
760 760 else:
761 761 timeLag = 15*10**-3
762 762 lag = numpy.ceil(timeLag/timeInterval)
763 763
764 764 for i in range(len(listMeteors)):
765 765
766 766 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
767 767 meteorAux = numpy.zeros(16)
768 768
769 769 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
770 770 mHeight = listMeteors[i][0]
771 771 mStart = listMeteors[i][1]
772 772 mPeak = listMeteors[i][2]
773 773 mEnd = listMeteors[i][3]
774 774
775 775 #get the volt data between the start and end times of the meteor
776 776 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
777 777 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
778 778
779 779 #3.6. Phase Difference estimation
780 780 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
781 781
782 782 #3.7. Phase difference removal & meteor start, peak and end times reestimated
783 783 #meteorVolts0.- all Channels, all Profiles
784 784 meteorVolts0 = volts[:,:,mHeight]
785 785 meteorThresh = noise[:,mHeight]*thresholdNoise
786 786 meteorNoise = noise[:,mHeight]
787 787 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
788 788 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
789 789
790 790 #Times reestimation
791 791 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
792 792 if mStart1.size > 0:
793 793 mStart1 = mStart1[-1] + 1
794 794
795 795 else:
796 796 mStart1 = mPeak
797 797
798 798 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
799 799 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
800 800 if mEndDecayTime1.size == 0:
801 801 mEndDecayTime1 = powerNet0.size
802 802 else:
803 803 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
804 804 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
805 805
806 806 #meteorVolts1.- all Channels, from start to end
807 807 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
808 808 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
809 809 if meteorVolts2.shape[1] == 0:
810 810 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
811 811 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
812 812 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
813 813 ##################### END PARAMETERS REESTIMATION #########################
814 814
815 815 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
816 816 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
817 817 if meteorVolts2.shape[1] > 0:
818 818 #Phase Difference re-estimation
819 819 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
820 820 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
821 821 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
822 822 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
823 823 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
824 824
825 825 #Phase Difference RMS
826 826 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
827 827 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
828 828 #Data from Meteor
829 829 mPeak1 = powerNet1.argmax() + mStart1
830 830 mPeakPower1 = powerNet1.max()
831 831 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
832 832 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
833 833 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
834 834 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
835 835 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
836 836 #Vectorize
837 837 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
838 838 meteorAux[7:11] = phaseDiffint[0:4]
839 839
840 840 #Rejection Criterions
841 841 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
842 842 meteorAux[-1] = 17
843 843 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
844 844 meteorAux[-1] = 1
845 845
846 846
847 847 else:
848 848 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
849 849 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
850 850 PowerSeries = 0
851 851
852 852 listMeteors1.append(meteorAux)
853 853 listPowerSeries.append(PowerSeries)
854 854 listVoltageSeries.append(meteorVolts1)
855 855
856 856 return listMeteors1, listPowerSeries, listVoltageSeries
857 857
858 858 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
859 859
860 860 threshError = 10
861 861 #Depending if it is 30 or 50 MHz
862 862 if frequency == 30e6:
863 863 timeLag = 45*10**-3
864 864 else:
865 865 timeLag = 15*10**-3
866 866 lag = numpy.ceil(timeLag/timeInterval)
867 867
868 868 listMeteors1 = []
869 869
870 870 for i in range(len(listMeteors)):
871 871 meteorPower = listPower[i]
872 872 meteorAux = listMeteors[i]
873 873
874 874 if meteorAux[-1] == 0:
875 875
876 876 try:
877 877 indmax = meteorPower.argmax()
878 878 indlag = indmax + lag
879 879
880 880 y = meteorPower[indlag:]
881 881 x = numpy.arange(0, y.size)*timeLag
882 882
883 883 #first guess
884 884 a = y[0]
885 885 tau = timeLag
886 886 #exponential fit
887 887 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
888 888 y1 = self.__exponential_function(x, *popt)
889 889 #error estimation
890 890 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
891 891
892 892 decayTime = popt[1]
893 893 riseTime = indmax*timeInterval
894 894 meteorAux[11:13] = [decayTime, error]
895 895
896 896 #Table items 7, 8 and 11
897 897 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
898 898 meteorAux[-1] = 7
899 899 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
900 900 meteorAux[-1] = 8
901 901 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
902 902 meteorAux[-1] = 11
903 903
904 904
905 905 except:
906 906 meteorAux[-1] = 11
907 907
908 908
909 909 listMeteors1.append(meteorAux)
910 910
911 911 return listMeteors1
912 912
913 913 #Exponential Function
914 914
915 915 def __exponential_function(self, x, a, tau):
916 916 y = a*numpy.exp(-x/tau)
917 917 return y
918 918
919 919 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
920 920
921 921 pairslist1 = list(pairslist)
922 922 pairslist1.append((0,1))
923 923 pairslist1.append((3,4))
924 924 numPairs = len(pairslist1)
925 925 #Time Lag
926 926 timeLag = 45*10**-3
927 927 c = 3e8
928 928 lag = numpy.ceil(timeLag/timeInterval)
929 929 freq = 30e6
930 930
931 931 listMeteors1 = []
932 932
933 933 for i in range(len(listMeteors)):
934 934 meteor = listMeteors[i]
935 935 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
936 936 if meteor[-1] == 0:
937 937 mStart = listMeteors[i][1]
938 938 mPeak = listMeteors[i][2]
939 939 mLag = mPeak - mStart + lag
940 940
941 941 #get the volt data between the start and end times of the meteor
942 942 meteorVolts = listVolts[i]
943 943 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
944 944
945 945 #Get CCF
946 946 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
947 947
948 948 #Method 2
949 949 slopes = numpy.zeros(numPairs)
950 950 time = numpy.array([-2,-1,1,2])*timeInterval
951 951 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
952 952
953 953 #Correct phases
954 954 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
955 955 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
956 956
957 957 if indDer[0].shape[0] > 0:
958 958 for i in range(indDer[0].shape[0]):
959 959 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
960 960 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
961 961
962 962 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
963 963 for j in range(numPairs):
964 964 fit = stats.linregress(time, angAllCCF[j,:])
965 965 slopes[j] = fit[0]
966 966
967 967 #Remove Outlier
968 968 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
969 969 # slopes = numpy.delete(slopes,indOut)
970 970 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
971 971 # slopes = numpy.delete(slopes,indOut)
972 972
973 973 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
974 974 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
975 975 meteorAux[-2] = radialError
976 976 meteorAux[-3] = radialVelocity
977 977
978 978 #Setting Error
979 979 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
980 980 if numpy.abs(radialVelocity) > 200:
981 981 meteorAux[-1] = 15
982 982 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
983 983 elif radialError > radialStdThresh:
984 984 meteorAux[-1] = 12
985 985
986 986 listMeteors1.append(meteorAux)
987 987 return listMeteors1
988 988
989 989 def __setNewArrays(self, listMeteors, date, heiRang):
990 990
991 991 #New arrays
992 992 arrayMeteors = numpy.array(listMeteors)
993 993 arrayParameters = numpy.zeros((len(listMeteors),10))
994 994
995 995 #Date inclusion
996 996 date = re.findall(r'\((.*?)\)', date)
997 997 date = date[0].split(',')
998 998 date = map(int, date)
999 999 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
1000 1000 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1001 1001
1002 1002 #Meteor array
1003 1003 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1004 1004 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1005 1005
1006 1006 #Parameters Array
1007 1007 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1008 1008 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1009 1009
1010 1010 return arrayMeteors, arrayParameters
1011 1011
1012 1012 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1013 1013
1014 1014 arrayAOA = numpy.zeros((phases.shape[0],3))
1015 1015 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1016 1016
1017 1017 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1018 1018 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1019 1019 arrayAOA[:,2] = cosDirError
1020 1020
1021 1021 azimuthAngle = arrayAOA[:,0]
1022 1022 zenithAngle = arrayAOA[:,1]
1023 1023
1024 1024 #Setting Error
1025 1025 #Number 3: AOA not fesible
1026 1026 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1027 1027 error[indInvalid] = 3
1028 1028 #Number 4: Large difference in AOAs obtained from different antenna baselines
1029 1029 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1030 1030 error[indInvalid] = 4
1031 1031 return arrayAOA, error
1032 1032
1033 1033 def __getDirectionCosines(self, arrayPhase, pairsList):
1034 1034
1035 1035 #Initializing some variables
1036 1036 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1037 1037 ang_aux = ang_aux.reshape(1,ang_aux.size)
1038 1038
1039 1039 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1040 1040 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1041 1041
1042 1042
1043 1043 for i in range(2):
1044 1044 #First Estimation
1045 1045 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1046 1046 #Dealias
1047 1047 indcsi = numpy.where(phi0_aux > numpy.pi)
1048 1048 phi0_aux[indcsi] -= 2*numpy.pi
1049 1049 indcsi = numpy.where(phi0_aux < -numpy.pi)
1050 1050 phi0_aux[indcsi] += 2*numpy.pi
1051 1051 #Direction Cosine 0
1052 1052 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1053 1053
1054 1054 #Most-Accurate Second Estimation
1055 1055 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1056 1056 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1057 1057 #Direction Cosine 1
1058 1058 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1059 1059
1060 1060 #Searching the correct Direction Cosine
1061 1061 cosdir0_aux = cosdir0[:,i]
1062 1062 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1063 1063 #Minimum Distance
1064 1064 cosDiff = (cosdir1 - cosdir0_aux)**2
1065 1065 indcos = cosDiff.argmin(axis = 1)
1066 1066 #Saving Value obtained
1067 1067 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1068 1068
1069 1069 return cosdir0, cosdir
1070 1070
1071 1071 def __calculateAOA(self, cosdir, azimuth):
1072 1072 cosdirX = cosdir[:,0]
1073 1073 cosdirY = cosdir[:,1]
1074 1074
1075 1075 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1076 1076 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1077 1077 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1078 1078
1079 1079 return angles
1080 1080
1081 1081 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1082 1082
1083 1083 Ramb = 375 #Ramb = c/(2*PRF)
1084 1084 Re = 6371 #Earth Radius
1085 1085 heights = numpy.zeros(Ranges.shape)
1086 1086
1087 1087 R_aux = numpy.array([0,1,2])*Ramb
1088 1088 R_aux = R_aux.reshape(1,R_aux.size)
1089 1089
1090 1090 Ranges = Ranges.reshape(Ranges.size,1)
1091 1091
1092 1092 Ri = Ranges + R_aux
1093 1093 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1094 1094
1095 1095 #Check if there is a height between 70 and 110 km
1096 1096 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1097 1097 ind_h = numpy.where(h_bool == 1)[0]
1098 1098
1099 1099 hCorr = hi[ind_h, :]
1100 1100 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1101 1101
1102 1102 hCorr = hi[ind_hCorr]
1103 1103 heights[ind_h] = hCorr
1104 1104
1105 1105 #Setting Error
1106 1106 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1107 1107 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1108 1108
1109 1109 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1110 1110 error[indInvalid2] = 14
1111 1111 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1112 1112 error[indInvalid1] = 13
1113 1113
1114 1114 return heights, error
1115 1115
1116 1116 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1117 1117
1118 1118 '''
1119 1119 Function GetMoments()
1120 1120
1121 1121 Input:
1122 1122 Output:
1123 1123 Variables modified:
1124 1124 '''
1125 1125 if path != None:
1126 1126 sys.path.append(path)
1127 1127 self.dataOut.library = importlib.import_module(file)
1128 1128
1129 1129 #To be inserted as a parameter
1130 1130 groupArray = numpy.array(groupList)
1131 1131 # groupArray = numpy.array([[0,1],[2,3]])
1132 1132 self.dataOut.groupList = groupArray
1133 1133
1134 1134 nGroups = groupArray.shape[0]
1135 1135 nChannels = self.dataIn.nChannels
1136 1136 nHeights=self.dataIn.heightList.size
1137 1137
1138 1138 #Parameters Array
1139 1139 self.dataOut.data_param = None
1140 1140
1141 1141 #Set constants
1142 1142 constants = self.dataOut.library.setConstants(self.dataIn)
1143 1143 self.dataOut.constants = constants
1144 1144 M = self.dataIn.normFactor
1145 1145 N = self.dataIn.nFFTPoints
1146 1146 ippSeconds = self.dataIn.ippSeconds
1147 1147 K = self.dataIn.nIncohInt
1148 1148 pairsArray = numpy.array(self.dataIn.pairsList)
1149 1149
1150 1150 #List of possible combinations
1151 1151 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1152 1152 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1153 1153
1154 1154 if getSNR:
1155 1155 listChannels = groupArray.reshape((groupArray.size))
1156 1156 listChannels.sort()
1157 1157 noise = self.dataIn.getNoise()
1158 self.dataOut.SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1158 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1159 1159
1160 1160 for i in range(nGroups):
1161 1161 coord = groupArray[i,:]
1162 1162
1163 1163 #Input data array
1164 1164 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1165 1165 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1166 1166
1167 1167 #Cross Spectra data array for Covariance Matrixes
1168 1168 ind = 0
1169 1169 for pairs in listComb:
1170 1170 pairsSel = numpy.array([coord[x],coord[y]])
1171 1171 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1172 1172 ind += 1
1173 1173 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1174 1174 dataCross = dataCross**2/K
1175 1175
1176 1176 for h in range(nHeights):
1177 1177 # print self.dataOut.heightList[h]
1178 1178
1179 1179 #Input
1180 1180 d = data[:,h]
1181 1181
1182 1182 #Covariance Matrix
1183 1183 D = numpy.diag(d**2/K)
1184 1184 ind = 0
1185 1185 for pairs in listComb:
1186 1186 #Coordinates in Covariance Matrix
1187 1187 x = pairs[0]
1188 1188 y = pairs[1]
1189 1189 #Channel Index
1190 1190 S12 = dataCross[ind,:,h]
1191 1191 D12 = numpy.diag(S12)
1192 1192 #Completing Covariance Matrix with Cross Spectras
1193 1193 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1194 1194 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1195 1195 ind += 1
1196 1196 Dinv=numpy.linalg.inv(D)
1197 1197 L=numpy.linalg.cholesky(Dinv)
1198 1198 LT=L.T
1199 1199
1200 1200 dp = numpy.dot(LT,d)
1201 1201
1202 1202 #Initial values
1203 1203 data_spc = self.dataIn.data_spc[coord,:,h]
1204 p0 = self.dataOut.library.initialValuesFunction(data_spc, constants)
1204 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants))
1205 1205
1206 try:
1206 1207 #Least Squares
1207 1208 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1208 1209 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1209 1210 #Chi square error
1210 1211 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1211 # error0 = 0
1212 1212 #Error with Jacobian
1213 1213 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1214 except:
1215 minp = p0*numpy.nan
1216 error0 = numpy.nan
1217 error1 = p0*numpy.nan
1218
1214 1219 #Save
1215 1220 if self.dataOut.data_param == None:
1216 self.dataOut.data_param = numpy.zeros((nGroups, minp.size, nHeights))*numpy.nan
1217 self.dataOut.error = numpy.zeros((nGroups, error1.size + 1, nHeights))*numpy.nan
1221 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1222 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1218 1223
1219 self.dataOut.error[i,:,h] = numpy.hstack((error0,error1))
1224 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1220 1225 self.dataOut.data_param[i,:,h] = minp
1221 1226 return
1222 1227
1223 1228
1224 1229 def __residFunction(self, p, dp, LT, constants):
1225 1230
1226 1231 fm = self.dataOut.library.modelFunction(p, constants)
1227 1232 fmp=numpy.dot(LT,fm)
1228 1233
1229 1234 return dp-fmp
1230 1235
1231 1236 def __getSNR(self, z, noise):
1232 1237
1233 1238 avg = numpy.average(z, axis=1)
1234 1239 SNR = (avg.T-noise)/noise
1235 1240 SNR = SNR.T
1236 1241 return SNR
1237 1242
1238 1243 def __chisq(p,chindex,hindex):
1239 1244 #similar to Resid but calculates CHI**2
1240 1245 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1241 1246 dp=numpy.dot(LT,d)
1242 1247 fmp=numpy.dot(LT,fm)
1243 1248 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1244 1249 return chisq
1245 1250
1246 1251
1247 1252
1248 1253 class WindProfiler(Operation):
1249 1254
1250 1255 __isConfig = False
1251 1256
1252 1257 __initime = None
1253 1258 __lastdatatime = None
1254 1259 __integrationtime = None
1255 1260
1256 1261 __buffer = None
1257 1262
1258 1263 __dataReady = False
1259 1264
1260 1265 __firstdata = None
1261 1266
1262 1267 n = None
1263 1268
1264 1269 def __init__(self):
1265 1270 Operation.__init__(self)
1266 1271
1267 1272 def __calculateCosDir(self, elev, azim):
1268 1273 zen = (90 - elev)*numpy.pi/180
1269 1274 azim = azim*numpy.pi/180
1270 1275 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1271 1276 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1272 1277
1273 1278 signX = numpy.sign(numpy.cos(azim))
1274 1279 signY = numpy.sign(numpy.sin(azim))
1275 1280
1276 1281 cosDirX = numpy.copysign(cosDirX, signX)
1277 1282 cosDirY = numpy.copysign(cosDirY, signY)
1278 1283 return cosDirX, cosDirY
1279 1284
1280 1285 def __calculateAngles(self, theta_x, theta_y, azimuth):
1281 1286
1282 1287 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1283 1288 zenith_arr = numpy.arccos(dir_cosw)
1284 1289 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1285 1290
1286 1291 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1287 1292 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1288 1293
1289 1294 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1290 1295
1291 1296 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1292 1297
1293 1298 #
1294 1299 if horOnly:
1295 1300 A = numpy.c_[dir_cosu,dir_cosv]
1296 1301 else:
1297 1302 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1298 1303 A = numpy.asmatrix(A)
1299 1304 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1300 1305
1301 1306 return A1
1302 1307
1303 1308 def __correctValues(self, heiRang, phi, velRadial, SNR):
1304 1309 listPhi = phi.tolist()
1305 1310 maxid = listPhi.index(max(listPhi))
1306 1311 minid = listPhi.index(min(listPhi))
1307 1312
1308 1313 rango = range(len(phi))
1309 1314 # rango = numpy.delete(rango,maxid)
1310 1315
1311 1316 heiRang1 = heiRang*math.cos(phi[maxid])
1312 1317 heiRangAux = heiRang*math.cos(phi[minid])
1313 1318 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1314 1319 heiRang1 = numpy.delete(heiRang1,indOut)
1315 1320
1316 1321 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1317 1322 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1318 1323
1319 1324 for i in rango:
1320 1325 x = heiRang*math.cos(phi[i])
1321 1326 y1 = velRadial[i,:]
1322 1327 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1323 1328
1324 1329 x1 = heiRang1
1325 1330 y11 = f1(x1)
1326 1331
1327 1332 y2 = SNR[i,:]
1328 1333 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1329 1334 y21 = f2(x1)
1330 1335
1331 1336 velRadial1[i,:] = y11
1332 1337 SNR1[i,:] = y21
1333 1338
1334 1339 return heiRang1, velRadial1, SNR1
1335 1340
1336 1341 def __calculateVelUVW(self, A, velRadial):
1337 1342
1338 1343 #Operacion Matricial
1339 1344 # velUVW = numpy.zeros((velRadial.shape[1],3))
1340 1345 # for ind in range(velRadial.shape[1]):
1341 1346 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1342 1347 # velUVW = velUVW.transpose()
1343 1348 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1344 1349 velUVW[:,:] = numpy.dot(A,velRadial)
1345 1350
1346 1351
1347 1352 return velUVW
1348 1353
1349 1354 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1350 1355 """
1351 1356 Function that implements Doppler Beam Swinging (DBS) technique.
1352 1357
1353 1358 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1354 1359 Direction correction (if necessary), Ranges and SNR
1355 1360
1356 1361 Output: Winds estimation (Zonal, Meridional and Vertical)
1357 1362
1358 1363 Parameters affected: Winds, height range, SNR
1359 1364 """
1360 1365 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1361 1366 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1362 1367 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1363 1368
1364 1369 #Calculo de Componentes de la velocidad con DBS
1365 1370 winds = self.__calculateVelUVW(A,velRadial1)
1366 1371
1367 1372 return winds, heiRang1, SNR1
1368 1373
1369 1374 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1370 1375
1371 1376 posx = numpy.asarray(posx)
1372 1377 posy = numpy.asarray(posy)
1373 1378
1374 1379 #Rotacion Inversa para alinear con el azimuth
1375 1380 if azimuth!= None:
1376 1381 azimuth = azimuth*math.pi/180
1377 1382 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1378 1383 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1379 1384 else:
1380 1385 posx1 = posx
1381 1386 posy1 = posy
1382 1387
1383 1388 #Calculo de Distancias
1384 1389 distx = numpy.zeros(pairsCrossCorr.size)
1385 1390 disty = numpy.zeros(pairsCrossCorr.size)
1386 1391 dist = numpy.zeros(pairsCrossCorr.size)
1387 1392 ang = numpy.zeros(pairsCrossCorr.size)
1388 1393
1389 1394 for i in range(pairsCrossCorr.size):
1390 1395 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1391 1396 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1392 1397 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1393 1398 ang[i] = numpy.arctan2(disty[i],distx[i])
1394 1399 #Calculo de Matrices
1395 1400 nPairs = len(pairs)
1396 1401 ang1 = numpy.zeros((nPairs, 2, 1))
1397 1402 dist1 = numpy.zeros((nPairs, 2, 1))
1398 1403
1399 1404 for j in range(nPairs):
1400 1405 dist1[j,0,0] = dist[pairs[j][0]]
1401 1406 dist1[j,1,0] = dist[pairs[j][1]]
1402 1407 ang1[j,0,0] = ang[pairs[j][0]]
1403 1408 ang1[j,1,0] = ang[pairs[j][1]]
1404 1409
1405 1410 return distx,disty, dist1,ang1
1406 1411
1407 1412 def __calculateVelVer(self, phase, lagTRange, _lambda):
1408 1413
1409 1414 Ts = lagTRange[1] - lagTRange[0]
1410 1415 velW = -_lambda*phase/(4*math.pi*Ts)
1411 1416
1412 1417 return velW
1413 1418
1414 1419 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1415 1420 nPairs = tau1.shape[0]
1416 1421 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1417 1422
1418 1423 angCos = numpy.cos(ang)
1419 1424 angSin = numpy.sin(ang)
1420 1425
1421 1426 vel0 = dist*tau1/(2*tau2**2)
1422 1427 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1423 1428 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1424 1429
1425 1430 ind = numpy.where(numpy.isinf(vel))
1426 1431 vel[ind] = numpy.nan
1427 1432
1428 1433 return vel
1429 1434
1430 1435 def __getPairsAutoCorr(self, pairsList, nChannels):
1431 1436
1432 1437 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1433 1438
1434 1439 for l in range(len(pairsList)):
1435 1440 firstChannel = pairsList[l][0]
1436 1441 secondChannel = pairsList[l][1]
1437 1442
1438 1443 #Obteniendo pares de Autocorrelacion
1439 1444 if firstChannel == secondChannel:
1440 1445 pairsAutoCorr[firstChannel] = int(l)
1441 1446
1442 1447 pairsAutoCorr = pairsAutoCorr.astype(int)
1443 1448
1444 1449 pairsCrossCorr = range(len(pairsList))
1445 1450 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1446 1451
1447 1452 return pairsAutoCorr, pairsCrossCorr
1448 1453
1449 1454 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1450 1455 """
1451 1456 Function that implements Spaced Antenna (SA) technique.
1452 1457
1453 1458 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1454 1459 Direction correction (if necessary), Ranges and SNR
1455 1460
1456 1461 Output: Winds estimation (Zonal, Meridional and Vertical)
1457 1462
1458 1463 Parameters affected: Winds
1459 1464 """
1460 1465 #Cross Correlation pairs obtained
1461 1466 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1462 1467 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1463 1468 pairsSelArray = numpy.array(pairsSelected)
1464 1469 pairs = []
1465 1470
1466 1471 #Wind estimation pairs obtained
1467 1472 for i in range(pairsSelArray.shape[0]/2):
1468 1473 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1469 1474 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1470 1475 pairs.append((ind1,ind2))
1471 1476
1472 1477 indtau = tau.shape[0]/2
1473 1478 tau1 = tau[:indtau,:]
1474 1479 tau2 = tau[indtau:-1,:]
1475 1480 tau1 = tau1[pairs,:]
1476 1481 tau2 = tau2[pairs,:]
1477 1482 phase1 = tau[-1,:]
1478 1483
1479 1484 #---------------------------------------------------------------------
1480 1485 #Metodo Directo
1481 1486 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1482 1487 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1483 1488 winds = stats.nanmean(winds, axis=0)
1484 1489 #---------------------------------------------------------------------
1485 1490 #Metodo General
1486 1491 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1487 1492 # #Calculo Coeficientes de Funcion de Correlacion
1488 1493 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1489 1494 # #Calculo de Velocidades
1490 1495 # winds = self.calculateVelUV(F,G,A,B,H)
1491 1496
1492 1497 #---------------------------------------------------------------------
1493 1498 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1494 1499 winds = correctFactor*winds
1495 1500 return winds
1496 1501
1497 1502 def __checkTime(self, currentTime, paramInterval, outputInterval):
1498 1503
1499 1504 dataTime = currentTime + paramInterval
1500 1505 deltaTime = dataTime - self.__initime
1501 1506
1502 1507 if deltaTime >= outputInterval or deltaTime < 0:
1503 1508 self.__dataReady = True
1504 1509 return
1505 1510
1506 1511 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1507 1512 '''
1508 1513 Function that implements winds estimation technique with detected meteors.
1509 1514
1510 1515 Input: Detected meteors, Minimum meteor quantity to wind estimation
1511 1516
1512 1517 Output: Winds estimation (Zonal and Meridional)
1513 1518
1514 1519 Parameters affected: Winds
1515 1520 '''
1516 1521 #Settings
1517 1522 nInt = (heightMax - heightMin)/2
1518 1523 winds = numpy.zeros((2,nInt))*numpy.nan
1519 1524
1520 1525 #Filter errors
1521 1526 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1522 1527 finalMeteor = arrayMeteor[error,:]
1523 1528
1524 1529 #Meteor Histogram
1525 1530 finalHeights = finalMeteor[:,3]
1526 1531 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1527 1532 nMeteorsPerI = hist[0]
1528 1533 heightPerI = hist[1]
1529 1534
1530 1535 #Sort of meteors
1531 1536 indSort = finalHeights.argsort()
1532 1537 finalMeteor2 = finalMeteor[indSort,:]
1533 1538
1534 1539 # Calculating winds
1535 1540 ind1 = 0
1536 1541 ind2 = 0
1537 1542
1538 1543 for i in range(nInt):
1539 1544 nMet = nMeteorsPerI[i]
1540 1545 ind1 = ind2
1541 1546 ind2 = ind1 + nMet
1542 1547
1543 1548 meteorAux = finalMeteor2[ind1:ind2,:]
1544 1549
1545 1550 if meteorAux.shape[0] >= meteorThresh:
1546 1551 vel = meteorAux[:, 7]
1547 1552 zen = meteorAux[:, 5]*numpy.pi/180
1548 1553 azim = meteorAux[:, 4]*numpy.pi/180
1549 1554
1550 1555 n = numpy.cos(zen)
1551 1556 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1552 1557 # l = m*numpy.tan(azim)
1553 1558 l = numpy.sin(zen)*numpy.sin(azim)
1554 1559 m = numpy.sin(zen)*numpy.cos(azim)
1555 1560
1556 1561 A = numpy.vstack((l, m)).transpose()
1557 1562 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1558 1563 windsAux = numpy.dot(A1, vel)
1559 1564
1560 1565 winds[0,i] = windsAux[0]
1561 1566 winds[1,i] = windsAux[1]
1562 1567
1563 1568 return winds, heightPerI[:-1]
1564 1569
1565 1570 def run(self, dataOut, technique, **kwargs):
1566 1571
1567 1572 param = dataOut.data_param
1568 1573 if dataOut.abscissaRange != None:
1569 1574 absc = dataOut.abscissaRange[:-1]
1570 1575 noise = dataOut.noise
1571 1576 heightRange = dataOut.getHeiRange()
1572 SNR = dataOut.SNR
1577 SNR = dataOut.data_SNR
1573 1578
1574 1579 if technique == 'DBS':
1575 1580
1576 1581 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1577 1582 theta_x = numpy.array(kwargs['dirCosx'])
1578 1583 theta_y = numpy.array(kwargs['dirCosy'])
1579 1584 else:
1580 1585 elev = numpy.array(kwargs['elevation'])
1581 1586 azim = numpy.array(kwargs['azimuth'])
1582 1587 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1583 1588 azimuth = kwargs['correctAzimuth']
1584 1589 if kwargs.has_key('horizontalOnly'):
1585 1590 horizontalOnly = kwargs['horizontalOnly']
1586 1591 else: horizontalOnly = False
1587 1592 if kwargs.has_key('correctFactor'):
1588 1593 correctFactor = kwargs['correctFactor']
1589 1594 else: correctFactor = 1
1590 1595 if kwargs.has_key('channelList'):
1591 1596 channelList = kwargs['channelList']
1592 1597 if len(channelList) == 2:
1593 1598 horizontalOnly = True
1594 1599 arrayChannel = numpy.array(channelList)
1595 1600 param = param[arrayChannel,:,:]
1596 1601 theta_x = theta_x[arrayChannel]
1597 1602 theta_y = theta_y[arrayChannel]
1598 1603
1599 1604 velRadial0 = param[:,1,:] #Radial velocity
1600 dataOut.data_output, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1605 dataOut.data_output, dataOut.heightRange, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1601 1606
1602 1607 elif technique == 'SA':
1603 1608
1604 1609 #Parameters
1605 1610 position_x = kwargs['positionX']
1606 1611 position_y = kwargs['positionY']
1607 1612 azimuth = kwargs['azimuth']
1608 1613
1609 1614 if kwargs.has_key('crosspairsList'):
1610 1615 pairs = kwargs['crosspairsList']
1611 1616 else:
1612 1617 pairs = None
1613 1618
1614 1619 if kwargs.has_key('correctFactor'):
1615 1620 correctFactor = kwargs['correctFactor']
1616 1621 else:
1617 1622 correctFactor = 1
1618 1623
1619 1624 tau = dataOut.data_param
1620 1625 _lambda = dataOut.C/dataOut.frequency
1621 1626 pairsList = dataOut.groupList
1622 1627 nChannels = dataOut.nChannels
1623 1628
1624 1629 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1625 1630 dataOut.initUtcTime = dataOut.ltctime
1626 1631 dataOut.outputInterval = dataOut.timeInterval
1627 1632
1628 1633 elif technique == 'Meteors':
1629 1634 dataOut.flagNoData = True
1630 1635 self.__dataReady = False
1631 1636
1632 1637 if kwargs.has_key('nHours'):
1633 1638 nHours = kwargs['nHours']
1634 1639 else:
1635 1640 nHours = 1
1636 1641
1637 1642 if kwargs.has_key('meteorsPerBin'):
1638 1643 meteorThresh = kwargs['meteorsPerBin']
1639 1644 else:
1640 1645 meteorThresh = 6
1641 1646
1642 1647 if kwargs.has_key('hmin'):
1643 1648 hmin = kwargs['hmin']
1644 1649 else: hmin = 70
1645 1650 if kwargs.has_key('hmax'):
1646 1651 hmax = kwargs['hmax']
1647 1652 else: hmax = 110
1648 1653
1649 1654 dataOut.outputInterval = nHours*3600
1650 1655
1651 1656 if self.__isConfig == False:
1652 1657 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1653 1658 #Get Initial LTC time
1654 1659 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1655 1660 self.__isConfig = True
1656 1661
1657 1662 if self.__buffer == None:
1658 1663 self.__buffer = dataOut.data_param
1659 1664 self.__firstdata = copy.copy(dataOut)
1660 1665
1661 1666 else:
1662 1667 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1663 1668
1664 1669 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1665 1670
1666 1671 if self.__dataReady:
1667 1672 dataOut.initUtcTime = self.__initime
1668 1673 self.__initime = self.__initime + dataOut.outputInterval #to erase time offset
1669 1674
1670 1675 dataOut.data_output, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1671 1676 dataOut.flagNoData = False
1672 1677 self.__buffer = None
1673 1678
1674 1679 return
1675 1680
1676 1681 class EWDriftsEstimation(Operation):
1677 1682
1678 1683
1679 1684 def __init__(self):
1680 1685 Operation.__init__(self)
1681 1686
1682 1687 def __correctValues(self, heiRang, phi, velRadial, SNR):
1683 1688 listPhi = phi.tolist()
1684 1689 maxid = listPhi.index(max(listPhi))
1685 1690 minid = listPhi.index(min(listPhi))
1686 1691
1687 1692 rango = range(len(phi))
1688 1693 # rango = numpy.delete(rango,maxid)
1689 1694
1690 1695 heiRang1 = heiRang*math.cos(phi[maxid])
1691 1696 heiRangAux = heiRang*math.cos(phi[minid])
1692 1697 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1693 1698 heiRang1 = numpy.delete(heiRang1,indOut)
1694 1699
1695 1700 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1696 1701 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1697 1702
1698 1703 for i in rango:
1699 1704 x = heiRang*math.cos(phi[i])
1700 1705 y1 = velRadial[i,:]
1701 1706 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1702 1707
1703 1708 x1 = heiRang1
1704 1709 y11 = f1(x1)
1705 1710
1706 1711 y2 = SNR[i,:]
1707 1712 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1708 1713 y21 = f2(x1)
1709 1714
1710 1715 velRadial1[i,:] = y11
1711 1716 SNR1[i,:] = y21
1712 1717
1713 1718 return heiRang1, velRadial1, SNR1
1714 1719
1715 1720 def run(self, dataOut, zenith, zenithCorrection):
1716 1721 heiRang = dataOut.heightList
1717 1722 velRadial = dataOut.data_param[:,3,:]
1718 SNR = dataOut.SNR
1723 SNR = dataOut.data_SNR
1719 1724
1720 1725 zenith = numpy.array(zenith)
1721 1726 zenith -= zenithCorrection
1722 1727 zenith *= numpy.pi/180
1723 1728
1724 1729 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1725 1730
1726 1731 alp = zenith[0]
1727 1732 bet = zenith[1]
1728 1733
1729 1734 w_w = velRadial1[0,:]
1730 1735 w_e = velRadial1[1,:]
1731 1736
1732 1737 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1733 1738 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1734 1739
1735 1740 winds = numpy.vstack((u,w))
1736 1741
1737 1742 dataOut.heightList = heiRang1
1738 1743 dataOut.data_output = winds
1739 dataOut.SNR = SNR1
1744 dataOut.data_SNR = SNR1
1740 1745
1741 1746 dataOut.initUtcTime = dataOut.ltctime
1742 1747 dataOut.outputInterval = dataOut.timeInterval
1743 1748 return
1744 1749
1745 1750
1746 1751
1747 1752
1748 1753
1749 1754 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,135 +1,140
1 1 # DIAS 19 Y 20 FEB 2014
2 2 # Comprobacion de Resultados DBS con SA
3 3
4 4 import os, sys
5 5
6 6 path = os.path.split(os.getcwd())[0]
7 7 sys.path.append(path)
8 8
9 9 from controller import *
10 10
11 11 desc = "DBS Experiment Test"
12 12 filename = "DBStest.xml"
13 13
14 14 controllerObj = Project()
15 15
16 16 controllerObj.setup(id = '191', name='test01', description=desc)
17 17
18 18 #Experimentos
19 19
20 20 path = '/host/Jicamarca/EW_Drifts/d2012248'
21 21 pathFigure = '/home/propietario/workspace/Graficos/drifts'
22 22
23 23
24 24 path = "/home/soporte/Data/drifts"
25 25 pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba'
26 pathFile = '/home/soporte/Data/drifts/HDF5'
26 27
27 xmin = 11.75
28 xmax = 14.75
28 xmin = 0
29 xmax = 24
29 30 #------------------------------------------------------------------------------------------------
30 31 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
31 32 path=path,
32 startDate='2012/01/01',
33 endDate='2012/12/31',
33 startDate='2012/09/06',
34 endDate='2012/09/06',
34 35 startTime='00:00:00',
35 36 endTime='23:59:59',
36 37 online=0,
37 38 walk=1)
38 39
39 40 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
40 41
41 42 #--------------------------------------------------------------------------------------------------
42 43
43 44 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
44 45
45 46 opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other')
46 47 opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist')
47 48
48 49 opObj11 = procUnitConfObj0.addOperation(name='filterByHeights')
49 50 opObj11.addParameter(name='window', value='3', format='int')
50 51
51 52 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
52 53 # opObj11.addParameter(name='code', value='1,-1', format='floatlist')
53 54 # opObj11.addParameter(name='nCode', value='2', format='int')
54 55 # opObj11.addParameter(name='nBaud', value='1', format='int')
55 56
56 57 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
57 58 procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int')
58 59 procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int')
59 60 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')#,(2,3)
60 61
61 62 opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
62 63 # # opObj11.addParameter(name='minHei', value='320.0', format='float')
63 64 # # opObj11.addParameter(name='maxHei', value='350.0', format='float')
64 65 opObj11.addParameter(name='minHei', value='200.0', format='float')
65 66 opObj11.addParameter(name='maxHei', value='600.0', format='float')
66 67
67 68 opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
68 69 opObj11.addParameter(name='channelList', value='0,1,2,3', format='intlist')
69 70
70 71 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
71 72 opObj11.addParameter(name='timeInterval', value='300.0', format='float')
72 73
73 74 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
74 75
75 76 # opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
76 77 # opObj14.addParameter(name='id', value='1', format='int')
77 78 # # opObj14.addParameter(name='wintitle', value='Con interf', format='str')
78 79 # opObj14.addParameter(name='save', value='1', format='bool')
79 80 # opObj14.addParameter(name='figpath', value=pathFigure, format='str')
80 81 # # opObj14.addParameter(name='zmin', value='5', format='int')
81 82 # opObj14.addParameter(name='zmax', value='30', format='int')
82 83 #
83 84 # opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
84 85 # opObj12.addParameter(name='id', value='2', format='int')
85 86 # opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
86 87 # opObj12.addParameter(name='save', value='1', format='bool')
87 88 # opObj12.addParameter(name='figpath', value = pathFigure, format='str')
88 89 # opObj12.addParameter(name='xmin', value=xmin, format='float')
89 90 # opObj12.addParameter(name='xmax', value=xmax, format='float')
90 91 # # opObj12.addParameter(name='zmin', value='5', format='int')
91 92 # opObj12.addParameter(name='zmax', value='30', format='int')
92 93
93 94 #--------------------------------------------------------------------------------------------------
94 95
95 96 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
96 97 opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting')
97 98 opObj20.addParameter(name='path', value='/home/soporte/workspace/RemoteSystemsTempFiles', format='str')
98 99 opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
99 100 opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
100 101
101 102 opObj11 = procUnitConfObj2.addOperation(name='SpectralFittingPlot', optype='other')
102 103 opObj11.addParameter(name='id', value='3', format='int')
103 104 opObj11.addParameter(name='wintitle', value='DopplerPlot', format='str')
104 105 opObj11.addParameter(name='cutHeight', value='350', format='int')
105 106 opObj11.addParameter(name='fit', value='1', format='int')#1--True/include fit
106 107 opObj11.addParameter(name='save', value='1', format='bool')
107 108 opObj11.addParameter(name='figpath', value = pathFigure, format='str')
108 109
110 opObj12 = procUnitConfObj2.addOperation(name='HDF5Writer', optype='other')
111 opObj12.addParameter(name='path', value=pathFile)
112 opObj12.addParameter(name='blocksPerFile', value='3', format='int')
113
109 114 opObj11 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
110 115 opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist')
111 116 opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float')
112 117
113 118 opObj23 = procUnitConfObj2.addOperation(name='EWDriftsPlot', optype='other')
114 119 opObj23.addParameter(name='id', value='4', format='int')
115 120 opObj23.addParameter(name='wintitle', value='EW Drifts', format='str')
116 121 opObj23.addParameter(name='save', value='1', format='bool')
117 122 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
118 123 opObj23.addParameter(name='zminZonal', value='-150', format='int')
119 124 opObj23.addParameter(name='zmaxZonal', value='150', format='int')
120 125 opObj23.addParameter(name='zminVertical', value='-30', format='float')
121 126 opObj23.addParameter(name='zmaxVertical', value='30', format='float')
122 127 opObj23.addParameter(name='SNR_1', value='1', format='bool')
123 128 opObj23.addParameter(name='SNRmax', value='5', format='int')
124 129 # opObj23.addParameter(name='SNRthresh', value='-50', format='float')
125 130 opObj23.addParameter(name='xmin', value=xmin, format='float')
126 131 opObj23.addParameter(name='xmax', value=xmax, format='float')
127 132 #--------------------------------------------------------------------------------------------------
128 133 print "Escribiendo el archivo XML"
129 134 controllerObj.writeXml(filename)
130 135 print "Leyendo el archivo XML"
131 136 controllerObj.readXml(filename)
132 137
133 138 controllerObj.createObjects()
134 139 controllerObj.connectObjects()
135 140 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,155 +1,169
1 1 # DIAS 19 Y 20 FEB 2014
2 2 # Comprobacion de Resultados DBS con SA
3 3
4 4 import os, sys
5 5
6 6 path = os.path.split(os.getcwd())[0]
7 7 sys.path.append(path)
8 8
9 9 from controller import *
10 10
11 11 desc = "DBS Experiment Test"
12 12 filename = "DBStest.xml"
13 13
14 14 controllerObj = Project()
15 15
16 16 controllerObj.setup(id = '191', name='test01', description=desc)
17 17
18 18 #Experimentos
19 19
20 20 #2014050 19 Feb 2014
21 21 # path = '/home/soporte/Documents/MST_Data/DBS/d2014050'
22 22 # pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/'
23 23 # xmin = '15.5'
24 24 # xmax = '23.99999999'
25 25 # startTime = '17:25:00'
26 26 # filehdf5 = "DBS_2014050.hdf5"
27 27
28 28 #2014051 20 Feb 2014
29 29 path = '/home/soporte/Data/MST/DBS/d2014051'
30 30 pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/'
31 31 xmin = '0'
32 32 xmax = '7.5'
33 33 startTime = '00:00:00'
34 34 filehdf5 = "DBS_2014051.hdf5"
35 35
36 36
37 37
38 38 #------------------------------------------------------------------------------------------------
39 39 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
40 40 path=path,
41 41 startDate='2014/01/31',
42 42 endDate='2014/03/31',
43 43 startTime=startTime,
44 44 endTime='23:59:59',
45 45 online=0,
46 46 delay=5,
47 47 walk=0)
48 48
49 49 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
50 50
51 51
52 52 #--------------------------------------------------------------------------------------------------
53 53
54 54 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
55 55
56 56 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
57 57
58 58 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
59 59 opObj11.addParameter(name='n', value='256', format='int')
60 60 # opObj11.addParameter(name='n', value='16', format='int')
61 61
62 62 opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
63 63 opObj11.addParameter(name='minIndex', value='10', format='float')
64 64 opObj11.addParameter(name='maxIndex', value='60', format='float')
65 65
66 66 #---------------------------------------------------------------------------------------------------
67 67
68 68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
69 69 procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int')
70 70 procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int')
71 71 # procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int')
72 72 procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
73 73
74 74 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
75 75 opObj11.addParameter(name='n', value='5', format='int')
76 76
77 77 opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
78 78 opObj14.addParameter(name='id', value='1', format='int')
79 79 opObj14.addParameter(name='wintitle', value='Con interf', format='str')
80 80 opObj14.addParameter(name='save', value='1', format='bool')
81 81 opObj14.addParameter(name='figpath', value=pathFigure, format='str')
82 82 opObj14.addParameter(name='zmin', value='5', format='int')
83 83 opObj14.addParameter(name='zmax', value='90', format='int')
84 84
85 85 opObj12 = procUnitConfObj1.addOperation(name='removeInterference')
86 86 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
87 87 opObj13.addParameter(name='mode', value='1', format='int')
88 88
89 89 opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
90 90 opObj12.addParameter(name='id', value='2', format='int')
91 91 opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
92 92 opObj12.addParameter(name='save', value='1', format='bool')
93 93 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
94 94 opObj12.addParameter(name='xmin', value=xmin, format='float')
95 95 opObj12.addParameter(name='xmax', value=xmax, format='float')
96 96 opObj12.addParameter(name='zmin', value='5', format='int')
97 97 opObj12.addParameter(name='zmax', value='90', format='int')
98 98
99 99 #--------------------------------------------------------------------------------------------------
100 100
101 101 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
102 102 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
103 103
104 104 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
105 105 opObj21.addParameter(name='id', value='3', format='int')
106 106 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
107 107 opObj21.addParameter(name='save', value='1', format='bool')
108 108 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
109 109 opObj21.addParameter(name='zmin', value='5', format='int')
110 110 opObj21.addParameter(name='zmax', value='90', format='int')
111 111
112 opObj21 = procUnitConfObj2.addOperation(name='RadialVelocityPlot', optype='other')
112 opObj21 = procUnitConfObj2.addOperation(name='ParametersPlot', optype='other')
113 113 opObj21.addParameter(name='id', value='5', format='int')
114 114 opObj21.addParameter(name='wintitle', value='Radial Velocity Plot', format='str')
115 115 opObj21.addParameter(name='save', value='1', format='bool')
116 116 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
117 117 opObj21.addParameter(name='SNRmin', value='-10', format='int')
118 118 opObj21.addParameter(name='SNRmax', value='60', format='int')
119 119 opObj21.addParameter(name='SNRthresh', value='0', format='float')
120 120 opObj21.addParameter(name='xmin', value=xmin, format='float')
121 121 opObj21.addParameter(name='xmax', value=xmax, format='float')
122 122
123 opObj21 = procUnitConfObj2.addOperation(name='ParametersPlot', optype='other')
124 opObj21.addParameter(name='id', value='6', format='int')
125 opObj21.addParameter(name='wintitle', value='Spectral width Plot', format='str')
126 opObj21.addParameter(name='save', value='1', format='bool')
127 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
128 opObj21.addParameter(name='SNRmin', value='-10', format='int')
129 opObj21.addParameter(name='SNRmax', value='60', format='int')
130 opObj21.addParameter(name='SNRthresh', value='0', format='float')
131 opObj21.addParameter(name='xmin', value=xmin, format='float')
132 opObj21.addParameter(name='xmax', value=xmax, format='float')
133 opObj21.addParameter(name='zmin', value=0, format='float')
134 opObj21.addParameter(name='paramIndex', value=2, format='int')
135 opObj21.addParameter(name='onlyPositive', value=1, format='bool')
136
123 137 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
124 138 opObj22.addParameter(name='technique', value='DBS', format='str')
125 139 opObj22.addParameter(name='correctAzimuth', value='51.06', format='float')
126 140 opObj22.addParameter(name='correctFactor', value='-1', format='float')
127 141 opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist')
128 142 opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist')
129 143 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
130 144 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
131 145
132 146 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
133 147 opObj23.addParameter(name='id', value='4', format='int')
134 148 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
135 149 opObj23.addParameter(name='save', value='1', format='bool')
136 150 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
137 151 opObj23.addParameter(name='zmin', value='-10', format='int')
138 152 opObj23.addParameter(name='zmax', value='10', format='int')
139 153 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
140 154 opObj23.addParameter(name='zmax_ver', value='80', format='float')
141 155 opObj23.addParameter(name='SNRmin', value='-10', format='int')
142 156 opObj23.addParameter(name='SNRmax', value='60', format='int')
143 157 opObj23.addParameter(name='SNRthresh', value='0', format='float')
144 158 opObj23.addParameter(name='xmin', value=xmin, format='float')
145 159 opObj23.addParameter(name='xmax', value=xmax, format='float')
146 160
147 161 #--------------------------------------------------------------------------------------------------
148 162 print "Escribiendo el archivo XML"
149 163 controllerObj.writeXml(filename)
150 164 print "Leyendo el archivo XML"
151 165 controllerObj.readXml(filename)
152 166
153 167 controllerObj.createObjects()
154 168 controllerObj.connectObjects()
155 169 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,139 +1,139
1 1 # DIAS 19 Y 20 FEB 2014
2 2 # Comprobacion de Resultados DBS con SA
3 3
4 4 import os, sys
5 5
6 6 path = os.path.split(os.getcwd())[0]
7 7 sys.path.append(path)
8 8
9 9 from controller import *
10 10
11 11 desc = "SA Experiment Test"
12 12 filename = "SA2014050.xml"
13 13
14 14 controllerObj = Project()
15 15
16 16 controllerObj.setup(id = '191', name='test01', description=desc)
17 17
18 18
19 19 #Experimentos
20 20
21 21 #2014050 19 Feb 2014
22 22 path = '/home/soporte/Data/MST/SA/d2014050'
23 pathFigure = '/home/soporte/workspace/Graficos/SA/new1/'
23 pathFigure = '/home/soporte/workspace/Graficos/SA/prueba1/'
24 24 xmin = '15.5'
25 25 xmax = '24'
26 26 startTime = '15:30:00'
27 27 filehdf5 = "SA_2014050.hdf5"
28 28
29 29 #2014051 20 Feb 2014
30 30 # path = '/home/soporte/Data/MST/SA/d2014051'
31 31 # pathFigure = '/home/soporte/workspace/Graficos/SA/new/'
32 32 # xmin = '0.0'
33 33 # xmax = '8.0'
34 34 # startTime = '00:00:00'
35 35 # filehdf5 = "SA_2014051.hdf5"
36 36
37 37 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
38 38 path=path,
39 39 startDate='2014/01/01',
40 40 endDate='2014/03/31',
41 41 startTime=startTime,
42 42 endTime='23:59:59',
43 43 online=0,
44 44 delay=5,
45 45 walk=0)
46 46
47 47 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
48 48
49 49
50 50 #--------------------------------------------------------------------------------------------------
51 51
52 52 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
53 53
54 54 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
55 55
56 56 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
57 57 opObj11.addParameter(name='n', value='600', format='int')
58 58 # opObj11.addParameter(name='n', value='10', format='int')
59 59
60 60 opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
61 61 opObj11.addParameter(name='minIndex', value='10', format='float')
62 62 opObj11.addParameter(name='maxIndex', value='60', format='float')
63 63 #---------------------------------------------------------------------------------------------------
64 64 procUnitConfObj1 = controllerObj.addProcUnit(datatype='CorrelationProc', inputId=procUnitConfObj0.getId())
65 65 # procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(1,0),(2,3)', format='pairsList')
66 66 procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(0,3),(0,2),(1,3),(1,2),(0,1),(2,3)', format='pairsList')
67 67 procUnitConfObj1.addParameter(name='fullT', value='1', format='bool')
68 68 procUnitConfObj1.addParameter(name='removeDC', value='1', format='bool')
69 69 #procUnitConfObj1.addParameter(name='lagT', value='0,1,2,3', format='intlist')
70 70
71 71 opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other')
72 72 opObj12.addParameter(name='id', value='1', format='int')
73 73 opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str')
74 74 opObj12.addParameter(name='save', value='1', format='bool')
75 75 opObj12.addParameter(name='zmin', value='0', format='int')
76 76 opObj12.addParameter(name='zmax', value='1', format='int')
77 77 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
78 78
79 79 opObj12 = procUnitConfObj1.addOperation(name='removeNoise')
80 80 opObj12.addParameter(name='mode', value='2', format='int')
81 81 opObj12 = procUnitConfObj1.addOperation(name='calculateNormFactor')
82 82
83 83 opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other')
84 84 opObj12.addParameter(name='id', value='2', format='int')
85 85 opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str')
86 86 opObj12.addParameter(name='save', value='1', format='bool')
87 87 opObj12.addParameter(name='zmin', value='0', format='int')
88 88 opObj12.addParameter(name='zmax', value='1', format='int')
89 89 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
90 90
91 91 #---------------------------------------------------------------------------------------------------
92 92 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
93 93 opObj20 = procUnitConfObj2.addOperation(name='GetLags')
94 94
95 95 opObj21 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
96 96 opObj21.addParameter(name='technique', value='SA', format='str')
97 97 # opObj21.addParameter(name='correctFactor', value='-1', format='float')
98 98 opObj21.addParameter(name='positionX', value='36,0,36,0', format='floatlist')
99 99 opObj21.addParameter(name='positionY', value='36,0,0,36', format='floatlist')
100 100 opObj21.addParameter(name='azimuth', value='51.06', format='float')
101 101 opObj21.addParameter(name='crosspairsList', value='(0,3),(0,2),(1,3),(1,2),(0,1),(2,3)', format='pairsList')#COrregir
102 102 #
103 103 opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
104 104 opObj22.addParameter(name='id', value='4', format='int')
105 105 opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
106 106 opObj22.addParameter(name='save', value='1', format='bool')
107 107 opObj22.addParameter(name='figpath', value = pathFigure, format='str')
108 108 opObj22.addParameter(name='zmin', value='-15', format='int')
109 109 opObj22.addParameter(name='zmax', value='15', format='int')
110 110 opObj22.addParameter(name='zmin_ver', value='-80', format='float')
111 111 opObj22.addParameter(name='zmax_ver', value='80', format='float')
112 112 opObj22.addParameter(name='SNRmin', value='-20', format='int')
113 113 opObj22.addParameter(name='SNRmax', value='40', format='int')
114 114 opObj22.addParameter(name='SNRthresh', value='-3.5', format='float')
115 115 opObj22.addParameter(name='xmin', value=xmin, format='float')
116 116 opObj22.addParameter(name='xmax', value=xmax, format='float')
117 117 # #-----------------------------------------------------------------------------------
118 118 #
119 119 # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId())
120 120 # procUnitConfObj2.addParameter(name='nFFTPoints', value='128', format='int')
121 121 # procUnitConfObj2.addParameter(name='nProfiles', value='128', format='int')
122 122 # procUnitConfObj2.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
123 123 #
124 124 # opObj22 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='other')
125 125 # opObj22.addParameter(name='id', value='5', format='int')
126 126 # opObj22.addParameter(name='wintitle', value='Spectra Plot', format='str')
127 127 # opObj22.addParameter(name='save', value='1', format='bool')
128 128 # opObj22.addParameter(name='figpath', value = pathFigure, format='str')
129 129
130 130 #-----------------------------------------------------------------------------------
131 131
132 132 print "Escribiendo el archivo XML"
133 133 controllerObj.writeXml(filename)
134 134 print "Leyendo el archivo XML"
135 135 controllerObj.readXml(filename)
136 136
137 137 controllerObj.createObjects()
138 138 controllerObj.connectObjects()
139 139 controllerObj.run() No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now