##// END OF EJS Templates
schain
RSS
merge con digital_rf y pull
José Chávez -
r1067:7abdb8296c68 merge
parent child
Show More
Side by Side Unified
Context file:
r1067:7abdb8296c68 Loading diff...:
Show file before | Show file after | Hide comments
@@ -9,56 +9,53 import datetime
9 import traceback
9 import traceback
10 import math
10 import math
11 import time
11 import time
12 from multiprocessing import Process, Queue, cpu_count
12 from multiprocessing import Process, cpu_count
13
14 import schainpy
15 import schainpy.admin
16 from schainpy.utils.log import logToFile
17
13
18 from xml.etree.ElementTree import ElementTree, Element, SubElement, tostring
14 from xml.etree.ElementTree import ElementTree, Element, SubElement, tostring
19 from xml.dom import minidom
15 from xml.dom import minidom
20
16
17 import schainpy
18 import schainpy.admin
21 from schainpy.model import *
19 from schainpy.model import *
22 from time import sleep
20 from schainpy.utils import log
23
21
22 DTYPES = {
23 'Voltage': '.r',
24 'Spectra': '.pdata'
25 }
24
26
27 def MPProject(project, n=cpu_count()):
28 '''
29 Project wrapper to run schain in n processes
30 '''
25
31
26 def prettify(elem):
32 rconf = project.getReadUnitObj()
27 """Return a pretty-printed XML string for the Element.
33 op = rconf.getOperationObj('run')
28 """
34 dt1 = op.getParameterValue('startDate')
29 rough_string = tostring(elem, 'utf-8')
35 dt2 = op.getParameterValue('endDate')
30 reparsed = minidom.parseString(rough_string)
31 return reparsed.toprettyxml(indent=" ")
32
33 def multiSchain(child, nProcess=cpu_count(), startDate=None, endDate=None, by_day=False):
34 skip = 0
35 cursor = 0
36 nFiles = None
37 processes = []
38 dt1 = datetime.datetime.strptime(startDate, '%Y/%m/%d')
39 dt2 = datetime.datetime.strptime(endDate, '%Y/%m/%d')
40 days = (dt2 - dt1).days
36 days = (dt2 - dt1).days
41
37
42 for day in range(days+1):
38 for day in range(days+1):
43 skip = 0
39 skip = 0
44 cursor = 0
40 cursor = 0
45 q = Queue()
46 processes = []
41 processes = []
47 dt = (dt1 + datetime.timedelta(day)).strftime('%Y/%m/%d')
42 dt = dt1 + datetime.timedelta(day)
48 firstProcess = Process(target=child, args=(cursor, skip, q, dt))
43 dt_str = dt.strftime('%Y/%m/%d')
49 firstProcess.start()
44 reader = JRODataReader()
50 if by_day:
45 paths, files = reader.searchFilesOffLine(path=rconf.path,
51 continue
46 startDate=dt,
52 nFiles = q.get()
47 endDate=dt,
48 ext=DTYPES[rconf.datatype])
49 nFiles = len(files)
53 if nFiles==0:
50 if nFiles == 0:
54 continue
51 continue
55 firstProcess.terminate()
52 skip = int(math.ceil(nFiles/n))
56 skip = int(math.ceil(nFiles/nProcess))
53 while nFiles > cursor*skip:
57 while True:
54 rconf.update(startDate=dt_str, endDate=dt_str, cursor=cursor,
58 processes.append(Process(target=child, args=(cursor, skip, q, dt)))
55 skip=skip)
59 processes[cursor].start()
56 p = project.clone()
60 if nFiles < cursor*skip:
57 p.start()
61 break
58 processes.append(p)
62 cursor += 1
59 cursor += 1
63
60
64 def beforeExit(exctype, value, trace):
61 def beforeExit(exctype, value, trace):
@@ -75,7 +72,6 def multiSchain(child, nProcess=cpu_count(), startDate=None, endDate=None, by_da
75
72
76 time.sleep(3)
73 time.sleep(3)
77
74
78
79 class ParameterConf():
75 class ParameterConf():
80
76
81 id = None
77 id = None
@@ -112,7 +108,7 class ParameterConf():
112 return self.__formated_value
108 return self.__formated_value
113
109
114 if value == '':
110 if value == '':
115 raise ValueError, "%s: This parameter value is empty" %self.name
111 raise ValueError, '%s: This parameter value is empty' %self.name
116
112
117 if format == 'list':
113 if format == 'list':
118 strList = value.split(',')
114 strList = value.split(',')
@@ -122,10 +118,10 class ParameterConf():
122 return self.__formated_value
118 return self.__formated_value
123
119
124 if format == 'intlist':
120 if format == 'intlist':
125 """
121 '''
126 Example:
122 Example:
127 value = (0,1,2)
123 value = (0,1,2)
128 """
124 '''
129
125
130 new_value = ast.literal_eval(value)
126 new_value = ast.literal_eval(value)
131
127
@@ -137,10 +133,10 class ParameterConf():
137 return self.__formated_value
133 return self.__formated_value
138
134
139 if format == 'floatlist':
135 if format == 'floatlist':
140 """
136 '''
141 Example:
137 Example:
142 value = (0.5, 1.4, 2.7)
138 value = (0.5, 1.4, 2.7)
143 """
139 '''
144
140
145 new_value = ast.literal_eval(value)
141 new_value = ast.literal_eval(value)
146
142
@@ -170,38 +166,38 class ParameterConf():
170 return self.__formated_value
166 return self.__formated_value
171
167
172 if format == 'pairslist':
168 if format == 'pairslist':
173 """
169 '''
174 Example:
170 Example:
175 value = (0,1),(1,2)
171 value = (0,1),(1,2)
176 """
172 '''
177
173
178 new_value = ast.literal_eval(value)
174 new_value = ast.literal_eval(value)
179
175
180 if type(new_value) not in (tuple, list):
176 if type(new_value) not in (tuple, list):
181 raise ValueError, "%s has to be a tuple or list of pairs" %value
177 raise ValueError, '%s has to be a tuple or list of pairs' %value
182
178
183 if type(new_value[0]) not in (tuple, list):
179 if type(new_value[0]) not in (tuple, list):
184 if len(new_value) != 2:
180 if len(new_value) != 2:
185 raise ValueError, "%s has to be a tuple or list of pairs" %value
181 raise ValueError, '%s has to be a tuple or list of pairs' %value
186 new_value = [new_value]
182 new_value = [new_value]
187
183
188 for thisPair in new_value:
184 for thisPair in new_value:
189 if len(thisPair) != 2:
185 if len(thisPair) != 2:
190 raise ValueError, "%s has to be a tuple or list of pairs" %value
186 raise ValueError, '%s has to be a tuple or list of pairs' %value
191
187
192 self.__formated_value = new_value
188 self.__formated_value = new_value
193
189
194 return self.__formated_value
190 return self.__formated_value
195
191
196 if format == 'multilist':
192 if format == 'multilist':
197 """
193 '''
198 Example:
194 Example:
199 value = (0,1,2),(3,4,5)
195 value = (0,1,2),(3,4,5)
200 """
196 '''
201 multiList = ast.literal_eval(value)
197 multiList = ast.literal_eval(value)
202
198
203 if type(multiList[0]) == int:
199 if type(multiList[0]) == int:
204 multiList = ast.literal_eval("(" + value + ")")
200 multiList = ast.literal_eval('(' + value + ')')
205
201
206 self.__formated_value = multiList
202 self.__formated_value = multiList
207
203
@@ -263,7 +259,7 class ParameterConf():
263
259
264 def printattr(self):
260 def printattr(self):
265
261
266 print "Parameter[%s]: name = %s, value = %s, format = %s" %(self.id, self.name, self.value, self.format)
262 print 'Parameter[%s]: name = %s, value = %s, format = %s' %(self.id, self.name, self.value, self.format)
267
263
268 class OperationConf():
264 class OperationConf():
269
265
@@ -372,6 +368,8 class OperationConf():
372
368
373 def addParameter(self, name, value, format='str'):
369 def addParameter(self, name, value, format='str'):
374
370
371 if value is None:
372 return None
375 id = self.__getNewId()
373 id = self.__getNewId()
376
374
377 parmConfObj = ParameterConf()
375 parmConfObj = ParameterConf()
@@ -431,7 +429,7 class OperationConf():
431
429
432 def printattr(self):
430 def printattr(self):
433
431
434 print "%s[%s]: name = %s, type = %s, priority = %s" %(self.ELEMENTNAME,
432 print '%s[%s]: name = %s, type = %s, priority = %s' %(self.ELEMENTNAME,
435 self.id,
433 self.id,
436 self.name,
434 self.name,
437 self.type,
435 self.type,
@@ -444,12 +442,11 class OperationConf():
444
442
445
443
446 if self.type == 'self':
444 if self.type == 'self':
447 raise ValueError, "This operation type cannot be created"
445 raise ValueError, 'This operation type cannot be created'
448
446
449 if self.type == 'plotter':
447 if self.type == 'plotter':
450 #Plotter(plotter_name)
451 if not plotter_queue:
448 if not plotter_queue:
452 raise ValueError, "plotter_queue is not defined. Use:\nmyProject = Project()\nmyProject.setPlotterQueue(plotter_queue)"
449 raise ValueError, 'plotter_queue is not defined. Use:\nmyProject = Project()\nmyProject.setPlotterQueue(plotter_queue)'
453
450
454 opObj = Plotter(self.name, plotter_queue)
451 opObj = Plotter(self.name, plotter_queue)
455
452
@@ -564,7 +561,7 class ProcUnitConf():
564
561
565 #Compatible with old signal chain version
562 #Compatible with old signal chain version
566 if datatype==None and name==None:
563 if datatype==None and name==None:
567 raise ValueError, "datatype or name should be defined"
564 raise ValueError, 'datatype or name should be defined'
568
565
569 if name==None:
566 if name==None:
570 if 'Proc' in datatype:
567 if 'Proc' in datatype:
@@ -595,7 +592,7 class ProcUnitConf():
595
592
596 def addParameter(self, **kwargs):
593 def addParameter(self, **kwargs):
597 '''
594 '''
598 Add parameters to "run" operation
595 Add parameters to 'run' operation
599 '''
596 '''
600 opObj = self.opConfObjList[0]
597 opObj = self.opConfObjList[0]
601
598
@@ -633,11 +630,11 class ProcUnitConf():
633 self.datatype = upElement.get('datatype')
630 self.datatype = upElement.get('datatype')
634 self.inputId = upElement.get('inputId')
631 self.inputId = upElement.get('inputId')
635
632
636 if self.ELEMENTNAME == "ReadUnit":
633 if self.ELEMENTNAME == 'ReadUnit':
637 self.datatype = self.datatype.replace("Reader", "")
634 self.datatype = self.datatype.replace('Reader', '')
638
635
639 if self.ELEMENTNAME == "ProcUnit":
636 if self.ELEMENTNAME == 'ProcUnit':
640 self.datatype = self.datatype.replace("Proc", "")
637 self.datatype = self.datatype.replace('Proc', '')
641
638
642 if self.inputId == 'None':
639 if self.inputId == 'None':
643 self.inputId = '0'
640 self.inputId = '0'
@@ -653,7 +650,7 class ProcUnitConf():
653
650
654 def printattr(self):
651 def printattr(self):
655
652
656 print "%s[%s]: name = %s, datatype = %s, inputId = %s" %(self.ELEMENTNAME,
653 print '%s[%s]: name = %s, datatype = %s, inputId = %s' %(self.ELEMENTNAME,
657 self.id,
654 self.id,
658 self.name,
655 self.name,
659 self.datatype,
656 self.datatype,
@@ -707,18 +704,10 class ProcUnitConf():
707
704
708 kwargs[parmConfObj.name] = parmConfObj.getValue()
705 kwargs[parmConfObj.name] = parmConfObj.getValue()
709
706
710 #ini = time.time()
711
712 #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id)
713 sts = self.procUnitObj.call(opType = opConfObj.type,
707 sts = self.procUnitObj.call(opType = opConfObj.type,
714 opName = opConfObj.name,
708 opName = opConfObj.name,
715 opId = opConfObj.id)
709 opId = opConfObj.id)
716
710
717 # total_time = time.time() - ini
718 #
719 # if total_time > 0.002:
720 # print "%s::%s took %f seconds" %(self.name, opConfObj.name, total_time)
721
722 is_ok = is_ok or sts
711 is_ok = is_ok or sts
723
712
724 return is_ok
713 return is_ok
@@ -762,11 +751,12 class ReadUnitConf(ProcUnitConf):
762
751
763 return self.ELEMENTNAME
752 return self.ELEMENTNAME
764
753
765 def setup(self, id, name, datatype, path='', startDate="", endDate="", startTime="",
754 def setup(self, id, name, datatype, path='', startDate='', endDate='',
766 endTime="", parentId=None, queue=None, server=None, **kwargs):
755 startTime='', endTime='', parentId=None, server=None, **kwargs):
756
767 #Compatible with old signal chain version
757 #Compatible with old signal chain version
768 if datatype==None and name==None:
758 if datatype==None and name==None:
769 raise ValueError, "datatype or name should be defined"
759 raise ValueError, 'datatype or name should be defined'
770
760
771 if name==None:
761 if name==None:
772 if 'Reader' in datatype:
762 if 'Reader' in datatype:
@@ -785,39 +775,28 class ReadUnitConf(ProcUnitConf):
785 self.endDate = endDate
775 self.endDate = endDate
786 self.startTime = startTime
776 self.startTime = startTime
787 self.endTime = endTime
777 self.endTime = endTime
788
789 self.inputId = '0'
778 self.inputId = '0'
790 self.parentId = parentId
779 self.parentId = parentId
791 self.queue = queue
792 self.server = server
780 self.server = server
793 self.addRunOperation(**kwargs)
781 self.addRunOperation(**kwargs)
794
782
795 def update(self, datatype, path, startDate, endDate, startTime, endTime, parentId=None, name=None, **kwargs):
783 def update(self, **kwargs):
796
797 #Compatible with old signal chain version
798 if datatype==None and name==None:
799 raise ValueError, "datatype or name should be defined"
800
784
801 if name==None:
785 if 'datatype' in kwargs:
786 datatype = kwargs.pop('datatype')
802 if 'Reader' in datatype:
787 if 'Reader' in datatype:
803 name = datatype
788 self.name = datatype
804 else:
789 else:
805 name = '%sReader' %(datatype)
790 self.name = '%sReader' %(datatype)
791 self.datatype = self.name.replace('Reader', '')
806
792
807 if datatype==None:
793 attrs = ('path', 'startDate', 'endDate', 'startTime', 'endTime', 'parentId')
808 datatype = name.replace('Reader','')
809
794
810 self.datatype = datatype
795 for attr in attrs:
811 self.name = name
796 if attr in kwargs:
812 self.path = path
797 setattr(self, attr, kwargs.pop(attr))
813 self.startDate = startDate
814 self.endDate = endDate
815 self.startTime = startTime
816 self.endTime = endTime
817
798
818 self.inputId = '0'
799 self.inputId = '0'
819 self.parentId = parentId
820
821 self.updateRunOperation(**kwargs)
800 self.updateRunOperation(**kwargs)
822
801
823 def removeOperations(self):
802 def removeOperations(self):
@@ -838,7 +817,7 class ReadUnitConf(ProcUnitConf):
838 opObj.addParameter(name='endDate' , value=self.endDate, format='date')
817 opObj.addParameter(name='endDate', value=self.endDate, format='date')
839 opObj.addParameter(name='startTime' , value=self.startTime, format='time')
818 opObj.addParameter(name='startTime', value=self.startTime, format='time')
840 opObj.addParameter(name='endTime' , value=self.endTime, format='time')
819 opObj.addParameter(name='endTime', value=self.endTime, format='time')
841 opObj.addParameter(name='queue' , value=self.queue, format='obj')
820
842 for key, value in kwargs.items():
821 for key, value in kwargs.items():
843 opObj.addParameter(name=key, value=value, format=type(value).__name__)
822 opObj.addParameter(name=key, value=value, format=type(value).__name__)
844 else:
823 else:
@@ -864,17 +843,6 class ReadUnitConf(ProcUnitConf):
864
843
865 return opObj
844 return opObj
866
845
867 # def makeXml(self, projectElement):
868 #
869 # procUnitElement = SubElement(projectElement, self.ELEMENTNAME)
870 # procUnitElement.set('id', str(self.id))
871 # procUnitElement.set('name', self.name)
872 # procUnitElement.set('datatype', self.datatype)
873 # procUnitElement.set('inputId', str(self.inputId))
874 #
875 # for opConfObj in self.opConfObjList:
876 # opConfObj.makeXml(procUnitElement)
877
878 def readXml(self, upElement):
846 def readXml(self, upElement):
879
847
880 self.id = upElement.get('id')
848 self.id = upElement.get('id')
@@ -882,8 +850,8 class ReadUnitConf(ProcUnitConf):
882 self.datatype = upElement.get('datatype')
850 self.datatype = upElement.get('datatype')
883 self.inputId = upElement.get('inputId')
851 self.inputId = upElement.get('inputId')
884
852
885 if self.ELEMENTNAME == "ReadUnit":
853 if self.ELEMENTNAME == 'ReadUnit':
886 self.datatype = self.datatype.replace("Reader", "")
854 self.datatype = self.datatype.replace('Reader', '')
887
855
888 if self.inputId == 'None':
856 if self.inputId == 'None':
889 self.inputId = '0'
857 self.inputId = '0'
@@ -905,8 +873,9 class ReadUnitConf(ProcUnitConf):
905 self.endTime = opConfObj.getParameterValue('endTime')
873 self.endTime = opConfObj.getParameterValue('endTime')
906
874
907 class Project(Process):
875 class Project(Process):
876
908 id = None
877 id = None
909 name = None
878 # name = None
910 description = None
879 description = None
911 filename = None
880 filename = None
912
881
@@ -916,13 +885,13 class Project(Process):
916
885
917 plotterQueue = None
886 plotterQueue = None
918
887
919 def __init__(self, plotter_queue=None, logfile=None):
888 def __init__(self, plotter_queue=None):
889
920 Process.__init__(self)
890 Process.__init__(self)
921 self.id = None
891 self.id = None
922 self.name = None
892 # self.name = None
923 self.description = None
893 self.description = None
924 if logfile is not None:
894
925 logToFile(logfile)
926 self.plotterQueue = plotter_queue
895 self.plotterQueue = plotter_queue
927
896
928 self.procUnitConfObjDict = {}
897 self.procUnitConfObjDict = {}
@@ -972,18 +941,28 class Project(Process):
972
941
973 self.procUnitConfObjDict = newProcUnitConfObjDict
942 self.procUnitConfObjDict = newProcUnitConfObjDict
974
943
975 def setup(self, id, name, description):
944 def setup(self, id, name='', description=''):
976
945
946 print
947 print '*'*60
948 print ' Starting SIGNAL CHAIN PROCESSING v%s ' % schainpy.__version__
949 print '*'*60
950 print
977 self.id = str(id)
951 self.id = str(id)
978 self.name = name
979 self.description = description
952 self.description = description
980
953
981 def update(self, name, description):
954 def update(self, name, description):
982
955
983 self.name = name
984 self.description = description
956 self.description = description
985
957
958 def clone(self):
959
960 p = Project()
961 p.procUnitConfObjDict = self.procUnitConfObjDict
962 return p
963
986 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
964 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
965
987 if id is None:
966 if id is None:
988 idReadUnit = self.__getNewId()
967 idReadUnit = self.__getNewId()
989 else:
968 else:
@@ -1021,7 +1000,7 class Project(Process):
1021 def getReadUnitObj(self):
1000 def getReadUnitObj(self):
1022
1001
1023 for obj in self.procUnitConfObjDict.values():
1002 for obj in self.procUnitConfObjDict.values():
1024 if obj.getElementName() == "ReadUnit":
1003 if obj.getElementName() == 'ReadUnit':
1025 return obj
1004 return obj
1026
1005
1027 return None
1006 return None
@@ -1066,20 +1045,20 class Project(Process):
1066 if self.filename:
1045 if self.filename:
1067 filename = self.filename
1046 filename = self.filename
1068 else:
1047 else:
1069 filename = "schain.xml"
1048 filename = 'schain.xml'
1070
1049
1071 if not filename:
1050 if not filename:
1072 print "filename has not been defined. Use setFilename(filename) for do it."
1051 print 'filename has not been defined. Use setFilename(filename) for do it.'
1073 return 0
1052 return 0
1074
1053
1075 abs_file = os.path.abspath(filename)
1054 abs_file = os.path.abspath(filename)
1076
1055
1077 if not os.access(os.path.dirname(abs_file), os.W_OK):
1056 if not os.access(os.path.dirname(abs_file), os.W_OK):
1078 print "No write permission on %s" %os.path.dirname(abs_file)
1057 print 'No write permission on %s' %os.path.dirname(abs_file)
1079 return 0
1058 return 0
1080
1059
1081 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
1060 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
1082 print "File %s already exists and it could not be overwriten" %abs_file
1061 print 'File %s already exists and it could not be overwriten' %abs_file
1083 return 0
1062 return 0
1084
1063
1085 self.makeXml()
1064 self.makeXml()
@@ -1093,13 +1072,13 class Project(Process):
1093 def readXml(self, filename = None):
1072 def readXml(self, filename = None):
1094
1073
1095 if not filename:
1074 if not filename:
1096 print "filename is not defined"
1075 print 'filename is not defined'
1097 return 0
1076 return 0
1098
1077
1099 abs_file = os.path.abspath(filename)
1078 abs_file = os.path.abspath(filename)
1100
1079
1101 if not os.path.isfile(abs_file):
1080 if not os.path.isfile(abs_file):
1102 print "%s file does not exist" %abs_file
1081 print '%s file does not exist' %abs_file
1103 return 0
1082 return 0
1104
1083
1105 self.projectElement = None
1084 self.projectElement = None
@@ -1108,7 +1087,7 class Project(Process):
1108 try:
1087 try:
1109 self.projectElement = ElementTree().parse(abs_file)
1088 self.projectElement = ElementTree().parse(abs_file)
1110 except:
1089 except:
1111 print "Error reading %s, verify file format" %filename
1090 print 'Error reading %s, verify file format' %filename
1112 return 0
1091 return 0
1113
1092
1114 self.project = self.projectElement.tag
1093 self.project = self.projectElement.tag
@@ -1145,7 +1124,7 class Project(Process):
1145
1124
1146 def printattr(self):
1125 def printattr(self):
1147
1126
1148 print "Project[%s]: name = %s, description = %s" %(self.id,
1127 print 'Project[%s]: name = %s, description = %s' %(self.id,
1149 self.name,
1128 self.name,
1150 self.description)
1129 self.description)
1151
1130
@@ -1179,7 +1158,7 class Project(Process):
1179
1158
1180 self.__connect(puObjIN, thisPUObj)
1159 self.__connect(puObjIN, thisPUObj)
1181
1160
1182 def __handleError(self, procUnitConfObj, send_email=True):
1161 def __handleError(self, procUnitConfObj, send_email=False):
1183
1162
1184 import socket
1163 import socket
1185
1164
@@ -1187,33 +1166,33 class Project(Process):
1187 sys.exc_info()[1],
1166 sys.exc_info()[1],
1188 sys.exc_info()[2])
1167 sys.exc_info()[2])
1189
1168
1190 print "***** Error occurred in %s *****" %(procUnitConfObj.name)
1169 print '***** Error occurred in %s *****' %(procUnitConfObj.name)
1191 print "***** %s" %err[-1]
1170 print '***** %s' %err[-1]
1192
1171
1193 message = "".join(err)
1172 message = ''.join(err)
1194
1173
1195 sys.stderr.write(message)
1174 sys.stderr.write(message)
1196
1175
1197 if not send_email:
1176 if not send_email:
1198 return
1177 return
1199
1178
1200 subject = "SChain v%s: Error running %s\n" %(schainpy.__version__, procUnitConfObj.name)
1179 subject = 'SChain v%s: Error running %s\n' %(schainpy.__version__, procUnitConfObj.name)
1201
1180
1202 subtitle = "%s: %s\n" %(procUnitConfObj.getElementName() ,procUnitConfObj.name)
1181 subtitle = '%s: %s\n' %(procUnitConfObj.getElementName() ,procUnitConfObj.name)
1203 subtitle += "Hostname: %s\n" %socket.gethostbyname(socket.gethostname())
1182 subtitle += 'Hostname: %s\n' %socket.gethostbyname(socket.gethostname())
1204 subtitle += "Working directory: %s\n" %os.path.abspath("./")
1183 subtitle += 'Working directory: %s\n' %os.path.abspath('./')
1205 subtitle += "Configuration file: %s\n" %self.filename
1184 subtitle += 'Configuration file: %s\n' %self.filename
1206 subtitle += "Time: %s\n" %str(datetime.datetime.now())
1185 subtitle += 'Time: %s\n' %str(datetime.datetime.now())
1207
1186
1208 readUnitConfObj = self.getReadUnitObj()
1187 readUnitConfObj = self.getReadUnitObj()
1209 if readUnitConfObj:
1188 if readUnitConfObj:
1210 subtitle += "\nInput parameters:\n"
1189 subtitle += '\nInput parameters:\n'
1211 subtitle += "[Data path = %s]\n" %readUnitConfObj.path
1190 subtitle += '[Data path = %s]\n' %readUnitConfObj.path
1212 subtitle += "[Data type = %s]\n" %readUnitConfObj.datatype
1191 subtitle += '[Data type = %s]\n' %readUnitConfObj.datatype
1213 subtitle += "[Start date = %s]\n" %readUnitConfObj.startDate
1192 subtitle += '[Start date = %s]\n' %readUnitConfObj.startDate
1214 subtitle += "[End date = %s]\n" %readUnitConfObj.endDate
1193 subtitle += '[End date = %s]\n' %readUnitConfObj.endDate
1215 subtitle += "[Start time = %s]\n" %readUnitConfObj.startTime
1194 subtitle += '[Start time = %s]\n' %readUnitConfObj.startTime
1216 subtitle += "[End time = %s]\n" %readUnitConfObj.endTime
1195 subtitle += '[End time = %s]\n' %readUnitConfObj.endTime
1217
1196
1218 adminObj = schainpy.admin.SchainNotify()
1197 adminObj = schainpy.admin.SchainNotify()
1219 adminObj.sendAlert(message=message,
1198 adminObj.sendAlert(message=message,
@@ -1228,15 +1207,15 class Project(Process):
1228 return 0
1207 return 0
1229
1208
1230 def runController(self):
1209 def runController(self):
1231 """
1210 '''
1232 returns 0 when this process has been stopped, 1 otherwise
1211 returns 0 when this process has been stopped, 1 otherwise
1233 """
1212 '''
1234
1213
1235 if self.isPaused():
1214 if self.isPaused():
1236 print "Process suspended"
1215 print 'Process suspended'
1237
1216
1238 while True:
1217 while True:
1239 sleep(0.1)
1218 time.sleep(0.1)
1240
1219
1241 if not self.isPaused():
1220 if not self.isPaused():
1242 break
1221 break
@@ -1244,10 +1223,10 class Project(Process):
1244 if self.isStopped():
1223 if self.isStopped():
1245 break
1224 break
1246
1225
1247 print "Process reinitialized"
1226 print 'Process reinitialized'
1248
1227
1249 if self.isStopped():
1228 if self.isStopped():
1250 print "Process stopped"
1229 print 'Process stopped'
1251 return 0
1230 return 0
1252
1231
1253 return 1
1232 return 1
@@ -1258,29 +1237,23 class Project(Process):
1258
1237
1259 def setPlotterQueue(self, plotter_queue):
1238 def setPlotterQueue(self, plotter_queue):
1260
1239
1261 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1240 raise NotImplementedError, 'Use schainpy.controller_api.ControllerThread instead Project class'
1262
1241
1263 def getPlotterQueue(self):
1242 def getPlotterQueue(self):
1264
1243
1265 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1244 raise NotImplementedError, 'Use schainpy.controller_api.ControllerThread instead Project class'
1266
1245
1267 def useExternalPlotter(self):
1246 def useExternalPlotter(self):
1268
1247
1269 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1248 raise NotImplementedError, 'Use schainpy.controller_api.ControllerThread instead Project class'
1270
1249
1250 def run(self):
1271
1251
1272 def run(self, filename=None):
1252 log.success('Starting {}'.format(self.name))
1273
1253
1274 # self.writeXml(filename)
1275 self.createObjects()
1254 self.createObjects()
1276 self.connectObjects()
1255 self.connectObjects()
1277
1256
1278 print
1279 print "*"*60
1280 print " Starting SIGNAL CHAIN PROCESSING v%s " %schainpy.__version__
1281 print "*"*60
1282 print
1283
1284 keyList = self.procUnitConfObjDict.keys()
1257 keyList = self.procUnitConfObjDict.keys()
1285 keyList.sort()
1258 keyList.sort()
1286
1259
@@ -1289,7 +1262,6 class Project(Process):
1289 is_ok = False
1262 is_ok = False
1290
1263
1291 for procKey in keyList:
1264 for procKey in keyList:
1292 # print "Running the '%s' process with %s" %(procUnitConfObj.name, procUnitConfObj.id)
1293
1265
1294 procUnitConfObj = self.procUnitConfObjDict[procKey]
1266 procUnitConfObj = self.procUnitConfObjDict[procKey]
1295
1267
@@ -1300,19 +1272,18 class Project(Process):
1300 is_ok = False
1272 is_ok = False
1301 break
1273 break
1302 except ValueError, e:
1274 except ValueError, e:
1303 sleep(0.5)
1275 time.sleep(0.5)
1304 self.__handleError(procUnitConfObj, send_email=True)
1276 self.__handleError(procUnitConfObj, send_email=True)
1305 is_ok = False
1277 is_ok = False
1306 break
1278 break
1307 except:
1279 except:
1308 sleep(0.5)
1280 time.sleep(0.5)
1309 self.__handleError(procUnitConfObj)
1281 self.__handleError(procUnitConfObj)
1310 is_ok = False
1282 is_ok = False
1311 break
1283 break
1312
1284
1313 #If every process unit finished so end process
1285 #If every process unit finished so end process
1314 if not(is_ok):
1286 if not(is_ok):
1315 # print "Every process unit have finished"
1316 break
1287 break
1317
1288
1318 if not self.runController():
1289 if not self.runController():
@@ -1322,3 +1293,5 class Project(Process):
1322 for procKey in keyList:
1293 for procKey in keyList:
1323 procUnitConfObj = self.procUnitConfObjDict[procKey]
1294 procUnitConfObj = self.procUnitConfObjDict[procKey]
1324 procUnitConfObj.close()
1295 procUnitConfObj.close()
1296
1297 log.success('{} finished'.format(self.name))
Show file before | Show file after | Hide comments
@@ -6218,3 +6218,6 class ShowMeConsole(QtCore.QObject):
6218 text = text[:-1]
6218 text = text[:-1]
6219
6219
6220 self.textWritten.emit(str(text))
6220 self.textWritten.emit(str(text))
6221
6222 def flush(self):
6223 pass
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
NO CONTENT: modified file
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1122 lines changed) Show them Hide them
@@ -1,32 +1,33
1
1
2 import os
2 import os
3 import zmq
4 import time
3 import time
5 import numpy
4 import glob
6 import datetime
5 import datetime
7 import numpy as np
6 from multiprocessing import Process
7
8 import zmq
9 import numpy
8 import matplotlib
10 import matplotlib
9 import glob
10 matplotlib.use('TkAgg')
11 import matplotlib.pyplot as plt
11 import matplotlib.pyplot as plt
12 from mpl_toolkits.axes_grid1 import make_axes_locatable
12 from mpl_toolkits.axes_grid1 import make_axes_locatable
13 from matplotlib.ticker import FuncFormatter, LinearLocator
13 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
14 from multiprocessing import Process
15
14
16 from schainpy.model.proc.jroproc_base import Operation
15 from schainpy.model.proc.jroproc_base import Operation
17
16 from schainpy.utils import log
18 plt.ion()
19
17
20 func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
18 func = lambda x, pos: ('%s') %(datetime.datetime.fromtimestamp(x).strftime('%H:%M'))
21 fromtimestamp = lambda x, mintime : (datetime.datetime.utcfromtimestamp(mintime).replace(hour=(x + 5), minute=0) - d1970).total_seconds()
22
23
19
24 d1970 = datetime.datetime(1970,1,1)
20 d1970 = datetime.datetime(1970, 1, 1)
25
21
22
26 class PlotData(Operation, Process):
23 class PlotData(Operation, Process):
24 '''
25 Base class for Schain plotting operations
26 '''
27
27
28 CODE = 'Figure'
28 CODE = 'Figure'
29 colormap = 'jro'
29 colormap = 'jro'
30 bgcolor = 'white'
30 CONFLATE = False
31 CONFLATE = False
31 __MAXNUMX = 80
32 __MAXNUMX = 80
32 __missing = 1E30
33 __missing = 1E30
@@ -37,54 +38,143 class PlotData(Operation, Process):
37 Process.__init__(self)
38 Process.__init__(self)
38 self.kwargs['code'] = self.CODE
39 self.kwargs['code'] = self.CODE
39 self.mp = False
40 self.mp = False
40 self.dataOut = None
41 self.data = None
41 self.isConfig = False
42 self.isConfig = False
42 self.figure = None
43 self.figures = []
43 self.axes = []
44 self.axes = []
45 self.cb_axes = []
44 self.localtime = kwargs.pop('localtime', True)
46 self.localtime = kwargs.pop('localtime', True)
45 self.show = kwargs.get('show', True)
47 self.show = kwargs.get('show', True)
46 self.save = kwargs.get('save', False)
48 self.save = kwargs.get('save', False)
47 self.colormap = kwargs.get('colormap', self.colormap)
49 self.colormap = kwargs.get('colormap', self.colormap)
48 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
50 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
49 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
51 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
50 self.showprofile = kwargs.get('showprofile', True)
52 self.colormaps = kwargs.get('colormaps', None)
51 self.title = kwargs.get('wintitle', '')
53 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
54 self.showprofile = kwargs.get('showprofile', False)
55 self.title = kwargs.get('wintitle', self.CODE.upper())
56 self.cb_label = kwargs.get('cb_label', None)
57 self.cb_labels = kwargs.get('cb_labels', None)
52 self.xaxis = kwargs.get('xaxis', 'frequency')
58 self.xaxis = kwargs.get('xaxis', 'frequency')
53 self.zmin = kwargs.get('zmin', None)
59 self.zmin = kwargs.get('zmin', None)
54 self.zmax = kwargs.get('zmax', None)
60 self.zmax = kwargs.get('zmax', None)
61 self.zlimits = kwargs.get('zlimits', None)
55 self.xmin = kwargs.get('xmin', None)
62 self.xmin = kwargs.get('xmin', None)
63 if self.xmin is not None:
64 self.xmin += 5
56 self.xmax = kwargs.get('xmax', None)
65 self.xmax = kwargs.get('xmax', None)
57 self.xrange = kwargs.get('xrange', 24)
66 self.xrange = kwargs.get('xrange', 24)
58 self.ymin = kwargs.get('ymin', None)
67 self.ymin = kwargs.get('ymin', None)
59 self.ymax = kwargs.get('ymax', None)
68 self.ymax = kwargs.get('ymax', None)
60 self.__MAXNUMY = kwargs.get('decimation', 5000)
69 self.xlabel = kwargs.get('xlabel', None)
61 self.throttle_value = 5
70 self.__MAXNUMY = kwargs.get('decimation', 100)
62 self.times = []
71 self.showSNR = kwargs.get('showSNR', False)
63 #self.interactive = self.kwargs['parent']
72 self.oneFigure = kwargs.get('oneFigure', True)
73 self.width = kwargs.get('width', None)
74 self.height = kwargs.get('height', None)
75 self.colorbar = kwargs.get('colorbar', True)
76 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
77 self.titles = ['' for __ in range(16)]
78
79 def __setup(self):
80 '''
81 Common setup for all figures, here figures and axes are created
82 '''
83
84 self.setup()
85
86 if self.width is None:
87 self.width = 8
88
89 self.figures = []
90 self.axes = []
91 self.cb_axes = []
92 self.pf_axes = []
93 self.cmaps = []
64
94
95 size = '15%' if self.ncols==1 else '30%'
96 pad = '4%' if self.ncols==1 else '8%'
97
98 if self.oneFigure:
99 if self.height is None:
100 self.height = 1.4*self.nrows + 1
101 fig = plt.figure(figsize=(self.width, self.height),
102 edgecolor='k',
103 facecolor='w')
104 self.figures.append(fig)
105 for n in range(self.nplots):
106 ax = fig.add_subplot(self.nrows, self.ncols, n+1)
107 ax.tick_params(labelsize=8)
108 ax.firsttime = True
109 self.axes.append(ax)
110 if self.showprofile:
111 cax = self.__add_axes(ax, size=size, pad=pad)
112 cax.tick_params(labelsize=8)
113 self.pf_axes.append(cax)
114 else:
115 if self.height is None:
116 self.height = 3
117 for n in range(self.nplots):
118 fig = plt.figure(figsize=(self.width, self.height),
119 edgecolor='k',
120 facecolor='w')
121 ax = fig.add_subplot(1, 1, 1)
122 ax.tick_params(labelsize=8)
123 ax.firsttime = True
124 self.figures.append(fig)
125 self.axes.append(ax)
126 if self.showprofile:
127 cax = self.__add_axes(ax, size=size, pad=pad)
128 cax.tick_params(labelsize=8)
129 self.pf_axes.append(cax)
130
131 for n in range(self.nrows):
132 if self.colormaps is not None:
133 cmap = plt.get_cmap(self.colormaps[n])
134 else:
135 cmap = plt.get_cmap(self.colormap)
136 cmap.set_bad(self.bgcolor, 1.)
137 self.cmaps.append(cmap)
138
139 def __add_axes(self, ax, size='30%', pad='8%'):
65 '''
140 '''
66 this new parameter is created to plot data from varius channels at different figures
141 Add new axes to the given figure
67 1. crear una lista de figuras donde se puedan plotear las figuras,
68 2. dar las opciones de configuracion a cada figura, estas opciones son iguales para ambas figuras
69 3. probar?
70 '''
142 '''
71 self.ind_plt_ch = kwargs.get('ind_plt_ch', False)
143 divider = make_axes_locatable(ax)
72 self.figurelist = None
144 nax = divider.new_horizontal(size=size, pad=pad)
145 ax.figure.add_axes(nax)
146 return nax
73
147
74
148
75 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
149 def setup(self):
150 '''
151 This method should be implemented in the child class, the following
152 attributes should be set:
76
153
154 self.nrows: number of rows
155 self.ncols: number of cols
156 self.nplots: number of plots (channels or pairs)
157 self.ylabel: label for Y axes
158 self.titles: list of axes title
159
160 '''
161 raise(NotImplementedError, 'Implement this method in child class')
162
163 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
164 '''
165 Create a masked array for missing data
166 '''
77 if x_buffer.shape[0] < 2:
167 if x_buffer.shape[0] < 2:
78 return x_buffer, y_buffer, z_buffer
168 return x_buffer, y_buffer, z_buffer
79
169
80 deltas = x_buffer[1:] - x_buffer[0:-1]
170 deltas = x_buffer[1:] - x_buffer[0:-1]
81 x_median = np.median(deltas)
171 x_median = numpy.median(deltas)
82
172
83 index = np.where(deltas > 5*x_median)
173 index = numpy.where(deltas > 5*x_median)
84
174
85 if len(index[0]) != 0:
175 if len(index[0]) != 0:
86 z_buffer[::, index[0], ::] = self.__missing
176 z_buffer[::, index[0], ::] = self.__missing
87 z_buffer = np.ma.masked_inside(z_buffer,
177 z_buffer = numpy.ma.masked_inside(z_buffer,
88 0.99*self.__missing,
178 0.99*self.__missing,
89 1.01*self.__missing)
179 1.01*self.__missing)
90
180
@@ -102,107 +192,114 class PlotData(Operation, Process):
102
192
103 return x, y, z
193 return x, y, z
104
194
195 def format(self):
105 '''
196 '''
106 JM:
197 Set min and max values, labels, ticks and titles
107 elimana las otras imagenes generadas debido a que lso workers no llegan en orden y le pueden
108 poner otro tiempo a la figura q no necesariamente es el ultimo.
109 Solo se realiza cuando termina la imagen.
110 Problemas:
111
112 File "/home/ci-81/workspace/schainv2.3/schainpy/model/graphics/jroplot_data.py", line 145, in __plot
113 for n, eachfigure in enumerate(self.figurelist):
114 TypeError: 'NoneType' object is not iterable
115
116 '''
198 '''
117 def deleteanotherfiles(self):
118 figurenames=[]
119 if self.figurelist != None:
120 for n, eachfigure in enumerate(self.figurelist):
121 #add specific name for each channel in channelList
122 ghostfigname = os.path.join(self.save, '{}_{}_{}'.format(self.titles[n].replace(' ',''),self.CODE,
123 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d')))
124 figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE,
125 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
126
199
127 for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures
200 if self.xmin is None:
128 if ghostfigure != figname:
201 xmin = self.min_time
129 os.remove(ghostfigure)
202 else:
130 print 'Removing GhostFigures:' , figname
203 if self.xaxis is 'time':
204 dt = datetime.datetime.fromtimestamp(self.min_time)
205 xmin = (datetime.datetime.combine(dt.date(),
206 datetime.time(int(self.xmin), 0, 0))-d1970).total_seconds()
131 else :
207 else:
132 '''Erasing ghost images for just on******************'''
208 xmin = self.xmin
133 ghostfigname = os.path.join(self.save, '{}_{}'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d')))
209
134 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
210 if self.xmax is None:
135 for ghostfigure in glob.glob(ghostfigname+'*'): #ghostfigure will adopt all posible names of figures
211 xmax = xmin+self.xrange*60*60
136 if ghostfigure != figname:
212 else:
137 os.remove(ghostfigure)
213 if self.xaxis is 'time':
138 print 'Removing GhostFigures:' , figname
214 dt = datetime.datetime.fromtimestamp(self.min_time)
215 xmax = (datetime.datetime.combine(dt.date(),
216 datetime.time(int(self.xmax), 0, 0))-d1970).total_seconds()
217 else:
218 xmax = self.xmax
219
220 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
221 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
222
223 ystep = 200 if ymax>= 800 else 100 if ymax>=400 else 50 if ymax>=200 else 20
224
225 for n, ax in enumerate(self.axes):
226 if ax.firsttime:
227 ax.set_facecolor(self.bgcolor)
228 ax.yaxis.set_major_locator(MultipleLocator(ystep))
229 if self.xaxis is 'time':
230 ax.xaxis.set_major_formatter(FuncFormatter(func))
231 ax.xaxis.set_major_locator(LinearLocator(9))
232 if self.xlabel is not None:
233 ax.set_xlabel(self.xlabel)
234 ax.set_ylabel(self.ylabel)
235 ax.firsttime = False
236 if self.showprofile:
237 self.pf_axes[n].set_ylim(ymin, ymax)
238 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
239 self.pf_axes[n].set_xlabel('dB')
240 self.pf_axes[n].grid(b=True, axis='x')
241 [tick.set_visible(False) for tick in self.pf_axes[n].get_yticklabels()]
242 if self.colorbar:
243 cb = plt.colorbar(ax.plt, ax=ax, pad=0.02)
244 cb.ax.tick_params(labelsize=8)
245 if self.cb_label:
246 cb.set_label(self.cb_label, size=8)
247 elif self.cb_labels:
248 cb.set_label(self.cb_labels[n], size=8)
249
250 ax.set_title('{} - {} UTC'.format(
251 self.titles[n],
252 datetime.datetime.fromtimestamp(self.max_time).strftime('%H:%M:%S')),
253 size=8)
254 ax.set_xlim(xmin, xmax)
255 ax.set_ylim(ymin, ymax)
256
139
257
140 def __plot(self):
258 def __plot(self):
259 '''
260 '''
261 log.success('Plotting', self.name)
141
262
142 print 'plotting...{}'.format(self.CODE)
143 if self.ind_plt_ch is False : #standard
144 if self.show:
145 self.figure.show()
146 self.plot()
263 self.plot()
147 plt.tight_layout()
264 self.format()
148 self.figure.canvas.manager.set_window_title('{} {} - {}'.format(self.title, self.CODE.upper(),
265
149 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
266 for n, fig in enumerate(self.figures):
150 else :
267 if self.nrows == 0 or self.nplots == 0:
151 print 'len(self.figurelist): ',len(self.figurelist)
268 log.warning('No data', self.name)
152 for n, eachfigure in enumerate(self.figurelist):
269 continue
153 if self.show:
270 if self.show:
154 eachfigure.show()
271 fig.show()
155
272
156 self.plot()
273 fig.tight_layout()
157 eachfigure.tight_layout() # ajuste de cada subplot
274 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
158 eachfigure.canvas.manager.set_window_title('{} {} - {}'.format(self.title[n], self.CODE.upper(),
159 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
275 datetime.datetime.fromtimestamp(self.max_time).strftime('%Y/%m/%d')))
276 # fig.canvas.draw()
160
277
161 # if self.save:
278 if self.save and self.data.ended:
162 # if self.ind_plt_ch is False : #standard
279 channels = range(self.nrows)
163 # figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
280 if self.oneFigure:
164 # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
281 label = ''
165 # print 'Saving figure: {}'.format(figname)
166 # self.figure.savefig(figname)
167 # else :
168 # for n, eachfigure in enumerate(self.figurelist):
169 # #add specific name for each channel in channelList
170 # figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n],self.CODE,
171 # datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
172 #
173 # print 'Saving figure: {}'.format(figname)
174 # eachfigure.savefig(figname)
175
176 if self.ind_plt_ch is False :
177 self.figure.canvas.draw()
178 else :
179 for eachfigure in self.figurelist:
180 eachfigure.canvas.draw()
181
182 if self.save:
183 if self.ind_plt_ch is False : #standard
184 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
185 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
186 print 'Saving figure: {}'.format(figname)
187 self.figure.savefig(figname)
188 else :
282 else:
189 for n, eachfigure in enumerate(self.figurelist):
283 label = '_{}'.format(channels[n])
190 #add specific name for each channel in channelList
284 figname = os.path.join(
191 figname = os.path.join(self.save, '{}_{}_{}.png'.format(self.titles[n].replace(' ',''),self.CODE,
285 self.save,
192 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')))
286 '{}{}_{}.png'.format(
193
287 self.CODE,
288 label,
289 datetime.datetime.fromtimestamp(self.saveTime).strftime('%y%m%d_%H%M%S')
290 )
291 )
194 print 'Saving figure: {}'.format(figname)
292 print 'Saving figure: {}'.format(figname)
195 eachfigure.savefig(figname)
293 fig.savefig(figname)
196
197
294
198 def plot(self):
295 def plot(self):
199
296 '''
200 print 'plotting...{}'.format(self.CODE.upper())
297 '''
201 return
298 raise(NotImplementedError, 'Implement this method in child class')
202
299
203 def run(self):
300 def run(self):
204
301
205 print '[Starting] {}'.format(self.name)
302 log.success('Starting', self.name)
206
303
207 context = zmq.Context()
304 context = zmq.Context()
208 receiver = context.socket(zmq.SUB)
305 receiver = context.socket(zmq.SUB)
@@ -214,150 +311,102 class PlotData(Operation, Process):
214 else:
311 else:
215 receiver.connect("ipc:///tmp/zmq.plots")
312 receiver.connect("ipc:///tmp/zmq.plots")
216
313
217 seconds_passed = 0
218
219 while True:
314 while True:
220 try:
315 try:
221 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)#flags=zmq.NOBLOCK
316 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
222 self.started = self.data['STARTED']
223 self.dataOut = self.data['dataOut']
224
317
225 if (len(self.times) < len(self.data['times']) and not self.started and self.data['ENDED']):
318 self.min_time = self.data.times[0]
226 continue
319 self.max_time = self.data.times[-1]
227
228 self.times = self.data['times']
229 self.times.sort()
230 self.throttle_value = self.data['throttle']
231 self.min_time = self.times[0]
232 self.max_time = self.times[-1]
233
320
234 if self.isConfig is False:
321 if self.isConfig is False:
235 print 'setting up'
322 self.__setup()
236 self.setup()
237 self.isConfig = True
323 self.isConfig = True
238 self.__plot()
239
324
240 if self.data['ENDED'] is True:
241 print '********GRAPHIC ENDED********'
242 self.ended = True
243 self.isConfig = False
244 self.__plot()
325 self.__plot()
245 self.deleteanotherfiles() #CLPDG
246 elif seconds_passed >= self.data['throttle']:
247 print 'passed', seconds_passed
248 self.__plot()
249 seconds_passed = 0
250
326
251 except zmq.Again as e:
327 except zmq.Again as e:
252 print 'Waiting for data...'
328 log.log('Waiting for data...')
253 plt.pause(2)
329 if self.data:
254 seconds_passed += 2
330 plt.pause(self.data.throttle)
331 else:
332 time.sleep(2)
255
333
256 def close(self):
334 def close(self):
257 if self.dataOut:
335 if self.data:
258 self.__plot()
336 self.__plot()
259
337
260
338
261 class PlotSpectraData(PlotData):
339 class PlotSpectraData(PlotData):
340 '''
341 Plot for Spectra data
342 '''
262
343
263 CODE = 'spc'
344 CODE = 'spc'
264 colormap = 'jro'
345 colormap = 'jro'
265 CONFLATE = False
266
346
267 def setup(self):
347 def setup(self):
268
348 self.nplots = len(self.data.channels)
269 ncolspan = 1
349 self.ncols = int(numpy.sqrt(self.nplots)+ 0.9)
270 colspan = 1
350 self.nrows = int((1.0*self.nplots/self.ncols) + 0.9)
271 self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
351 self.width = 3.4*self.ncols
272 self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
352 self.height = 3*self.nrows
273 self.width = 3.6*self.ncols
353 self.cb_label = 'dB'
274 self.height = 3.2*self.nrows
275 if self.showprofile:
354 if self.showprofile:
276 ncolspan = 3
355 self.width += 0.8*self.ncols
277 colspan = 2
278 self.width += 1.2*self.ncols
279
356
280 self.ylabel = 'Range [Km]'
357 self.ylabel = 'Range [Km]'
281 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
282
283 if self.figure is None:
284 self.figure = plt.figure(figsize=(self.width, self.height),
285 edgecolor='k',
286 facecolor='w')
287 else:
288 self.figure.clf()
289
290 n = 0
291 for y in range(self.nrows):
292 for x in range(self.ncols):
293 if n >= self.dataOut.nChannels:
294 break
295 ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
296 if self.showprofile:
297 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
298
299 ax.firsttime = True
300 self.axes.append(ax)
301 n += 1
302
358
303 def plot(self):
359 def plot(self):
304
305 if self.xaxis == "frequency":
360 if self.xaxis == "frequency":
306 x = self.dataOut.getFreqRange(1)/1000.
361 x = self.data.xrange[0]
307 xlabel = "Frequency (kHz)"
362 self.xlabel = "Frequency (kHz)"
308 elif self.xaxis == "time":
363 elif self.xaxis == "time":
309 x = self.dataOut.getAcfRange(1)
364 x = self.data.xrange[1]
310 xlabel = "Time (ms)"
365 self.xlabel = "Time (ms)"
311 else:
366 else:
312 x = self.dataOut.getVelRange(1)
367 x = self.data.xrange[2]
313 xlabel = "Velocity (m/s)"
368 self.xlabel = "Velocity (m/s)"
369
370 if self.CODE == 'spc_mean':
371 x = self.data.xrange[2]
372 self.xlabel = "Velocity (m/s)"
314
373
315 y = self.dataOut.getHeiRange()
374 self.titles = []
316 z = self.data[self.CODE]
375
376 y = self.data.heights
377 self.y = y
378 z = self.data['spc']
317
379
318 for n, ax in enumerate(self.axes):
380 for n, ax in enumerate(self.axes):
381 noise = self.data['noise'][n][-1]
382 if self.CODE == 'spc_mean':
383 mean = self.data['mean'][n][-1]
319 if ax.firsttime:
384 if ax.firsttime:
320 self.xmax = self.xmax if self.xmax else np.nanmax(x)
385 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
321 self.xmin = self.xmin if self.xmin else -self.xmax
386 self.xmin = self.xmin if self.xmin else -self.xmax
322 self.ymin = self.ymin if self.ymin else np.nanmin(y)
387 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
323 self.ymax = self.ymax if self.ymax else np.nanmax(y)
388 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
324 self.zmin = self.zmin if self.zmin else np.nanmin(z)
389 ax.plt = ax.pcolormesh(x, y, z[n].T,
325 self.zmax = self.zmax if self.zmax else np.nanmax(z)
326 ax.plot = ax.pcolormesh(x, y, z[n].T,
327 vmin=self.zmin,
390 vmin=self.zmin,
328 vmax=self.zmax,
391 vmax=self.zmax,
329 cmap=plt.get_cmap(self.colormap)
392 cmap=plt.get_cmap(self.colormap)
330 )
393 )
331 divider = make_axes_locatable(ax)
332 cax = divider.new_horizontal(size='3%', pad=0.05)
333 self.figure.add_axes(cax)
334 plt.colorbar(ax.plot, cax)
335
336 ax.set_xlim(self.xmin, self.xmax)
337 ax.set_ylim(self.ymin, self.ymax)
338
339 ax.set_ylabel(self.ylabel)
340 ax.set_xlabel(xlabel)
341
342 ax.firsttime = False
343
394
344 if self.showprofile:
395 if self.showprofile:
345 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
396 ax.plt_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], y)[0]
346 ax.ax_profile.set_xlim(self.zmin, self.zmax)
397 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
347 ax.ax_profile.set_ylim(self.ymin, self.ymax)
398 color="k", linestyle="dashed", lw=1)[0]
348 ax.ax_profile.set_xlabel('dB')
399 if self.CODE == 'spc_mean':
349 ax.ax_profile.grid(b=True, axis='x')
400 ax.plt_mean = ax.plot(mean, y, color='k')[0]
350 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
351 color="k", linestyle="dashed", lw=2)[0]
352 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
353 else:
401 else:
354 ax.plot.set_array(z[n].T.ravel())
402 ax.plt.set_array(z[n].T.ravel())
355 if self.showprofile:
403 if self.showprofile:
356 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
404 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
357 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
405 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
406 if self.CODE == 'spc_mean':
407 ax.plt_mean.set_data(mean, y)
358
408
359 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
409 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
360 size=8)
361 self.saveTime = self.max_time
410 self.saveTime = self.max_time
362
411
363
412
@@ -368,544 +417,244 class PlotCrossSpectraData(PlotData):
368 zmax_coh = None
417 zmax_coh = None
369 zmin_phase = None
418 zmin_phase = None
370 zmax_phase = None
419 zmax_phase = None
371 CONFLATE = False
372
420
373 def setup(self):
421 def setup(self):
374
422
375 ncolspan = 1
423 self.ncols = 4
376 colspan = 1
424 self.nrows = len(self.data.pairs)
377 self.ncols = 2
425 self.nplots = self.nrows*4
378 self.nrows = self.dataOut.nPairs
426 self.width = 3.4*self.ncols
379 self.width = 3.6*self.ncols
427 self.height = 3*self.nrows
380 self.height = 3.2*self.nrows
381
382 self.ylabel = 'Range [Km]'
428 self.ylabel = 'Range [Km]'
383 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
429 self.showprofile = False
384
385 if self.figure is None:
386 self.figure = plt.figure(figsize=(self.width, self.height),
387 edgecolor='k',
388 facecolor='w')
389 else:
390 self.figure.clf()
391
392 for y in range(self.nrows):
393 for x in range(self.ncols):
394 ax = plt.subplot2grid((self.nrows, self.ncols), (y, x), 1, 1)
395 ax.firsttime = True
396 self.axes.append(ax)
397
430
398 def plot(self):
431 def plot(self):
399
432
400 if self.xaxis == "frequency":
433 if self.xaxis == "frequency":
401 x = self.dataOut.getFreqRange(1)/1000.
434 x = self.data.xrange[0]
402 xlabel = "Frequency (kHz)"
435 self.xlabel = "Frequency (kHz)"
403 elif self.xaxis == "time":
436 elif self.xaxis == "time":
404 x = self.dataOut.getAcfRange(1)
437 x = self.data.xrange[1]
405 xlabel = "Time (ms)"
438 self.xlabel = "Time (ms)"
406 else:
439 else:
407 x = self.dataOut.getVelRange(1)
440 x = self.data.xrange[2]
408 xlabel = "Velocity (m/s)"
441 self.xlabel = "Velocity (m/s)"
442
443 self.titles = []
409
444
410 y = self.dataOut.getHeiRange()
445 y = self.data.heights
411 z_coh = self.data['cspc_coh']
446 self.y = y
412 z_phase = self.data['cspc_phase']
447 spc = self.data['spc']
448 cspc = self.data['cspc']
413
449
414 for n in range(self.nrows):
450 for n in range(self.nrows):
415 ax = self.axes[2*n]
451 noise = self.data['noise'][n][-1]
416 ax1 = self.axes[2*n+1]
452 pair = self.data.pairs[n]
453 ax = self.axes[4*n]
454 ax3 = self.axes[4*n+3]
417 if ax.firsttime:
455 if ax.firsttime:
418 self.xmax = self.xmax if self.xmax else np.nanmax(x)
456 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
419 self.xmin = self.xmin if self.xmin else -self.xmax
457 self.xmin = self.xmin if self.xmin else -self.xmax
420 self.ymin = self.ymin if self.ymin else np.nanmin(y)
458 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
421 self.ymax = self.ymax if self.ymax else np.nanmax(y)
459 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
422 self.zmin_coh = self.zmin_coh if self.zmin_coh else 0.0
460 ax.plt = ax.pcolormesh(x, y, spc[pair[0]].T,
423 self.zmax_coh = self.zmax_coh if self.zmax_coh else 1.0
461 vmin=self.zmin,
424 self.zmin_phase = self.zmin_phase if self.zmin_phase else -180
462 vmax=self.zmax,
425 self.zmax_phase = self.zmax_phase if self.zmax_phase else 180
463 cmap=plt.get_cmap(self.colormap)
426
427 ax.plot = ax.pcolormesh(x, y, z_coh[n].T,
428 vmin=self.zmin_coh,
429 vmax=self.zmax_coh,
430 cmap=plt.get_cmap(self.colormap_coh)
431 )
432 divider = make_axes_locatable(ax)
433 cax = divider.new_horizontal(size='3%', pad=0.05)
434 self.figure.add_axes(cax)
435 plt.colorbar(ax.plot, cax)
436
437 ax.set_xlim(self.xmin, self.xmax)
438 ax.set_ylim(self.ymin, self.ymax)
439
440 ax.set_ylabel(self.ylabel)
441 ax.set_xlabel(xlabel)
442 ax.firsttime = False
443
444 ax1.plot = ax1.pcolormesh(x, y, z_phase[n].T,
445 vmin=self.zmin_phase,
446 vmax=self.zmax_phase,
447 cmap=plt.get_cmap(self.colormap_phase)
448 )
464 )
449 divider = make_axes_locatable(ax1)
450 cax = divider.new_horizontal(size='3%', pad=0.05)
451 self.figure.add_axes(cax)
452 plt.colorbar(ax1.plot, cax)
453
454 ax1.set_xlim(self.xmin, self.xmax)
455 ax1.set_ylim(self.ymin, self.ymax)
456
457 ax1.set_ylabel(self.ylabel)
458 ax1.set_xlabel(xlabel)
459 ax1.firsttime = False
460 else:
465 else:
461 ax.plot.set_array(z_coh[n].T.ravel())
466 ax.plt.set_array(spc[pair[0]].T.ravel())
462 ax1.plot.set_array(z_phase[n].T.ravel())
467 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
463
464 ax.set_title('Coherence Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
465 ax1.set_title('Phase Ch{} * Ch{}'.format(self.dataOut.pairsList[n][0], self.dataOut.pairsList[n][1]), size=8)
466 self.saveTime = self.max_time
467
468
469 class PlotSpectraMeanData(PlotSpectraData):
470
471 CODE = 'spc_mean'
472 colormap = 'jet'
473
474 def plot(self):
475
476 if self.xaxis == "frequency":
477 x = self.dataOut.getFreqRange(1)/1000.
478 xlabel = "Frequency (kHz)"
479 elif self.xaxis == "time":
480 x = self.dataOut.getAcfRange(1)
481 xlabel = "Time (ms)"
482 else:
483 x = self.dataOut.getVelRange(1)
484 xlabel = "Velocity (m/s)"
485
486 y = self.dataOut.getHeiRange()
487 z = self.data['spc']
488 mean = self.data['mean'][self.max_time]
489
490 for n, ax in enumerate(self.axes):
491
468
469 ax = self.axes[4*n+1]
492 if ax.firsttime:
470 if ax.firsttime:
493 self.xmax = self.xmax if self.xmax else np.nanmax(x)
471 ax.plt = ax.pcolormesh(x, y, spc[pair[1]].T,
494 self.xmin = self.xmin if self.xmin else -self.xmax
495 self.ymin = self.ymin if self.ymin else np.nanmin(y)
496 self.ymax = self.ymax if self.ymax else np.nanmax(y)
497 self.zmin = self.zmin if self.zmin else np.nanmin(z)
498 self.zmax = self.zmax if self.zmax else np.nanmax(z)
499 ax.plt = ax.pcolormesh(x, y, z[n].T,
500 vmin=self.zmin,
472 vmin=self.zmin,
501 vmax=self.zmax,
473 vmax=self.zmax,
502 cmap=plt.get_cmap(self.colormap)
474 cmap=plt.get_cmap(self.colormap)
503 )
475 )
504 ax.plt_dop = ax.plot(mean[n], y,
476 else:
505 color='k')[0]
477 ax.plt.set_array(spc[pair[1]].T.ravel())
506
478 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
507 divider = make_axes_locatable(ax)
508 cax = divider.new_horizontal(size='3%', pad=0.05)
509 self.figure.add_axes(cax)
510 plt.colorbar(ax.plt, cax)
511
512 ax.set_xlim(self.xmin, self.xmax)
513 ax.set_ylim(self.ymin, self.ymax)
514
479
515 ax.set_ylabel(self.ylabel)
480 out = cspc[n]/numpy.sqrt(spc[pair[0]]*spc[pair[1]])
516 ax.set_xlabel(xlabel)
481 coh = numpy.abs(out)
482 phase = numpy.arctan2(out.imag, out.real)*180/numpy.pi
517
483
518 ax.firsttime = False
484 ax = self.axes[4*n+2]
485 if ax.firsttime:
486 ax.plt = ax.pcolormesh(x, y, coh.T,
487 vmin=0,
488 vmax=1,
489 cmap=plt.get_cmap(self.colormap_coh)
490 )
491 else:
492 ax.plt.set_array(coh.T.ravel())
493 self.titles.append('Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
519
494
520 if self.showprofile:
495 ax = self.axes[4*n+3]
521 ax.plt_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
496 if ax.firsttime:
522 ax.ax_profile.set_xlim(self.zmin, self.zmax)
497 ax.plt = ax.pcolormesh(x, y, phase.T,
523 ax.ax_profile.set_ylim(self.ymin, self.ymax)
498 vmin=-180,
524 ax.ax_profile.set_xlabel('dB')
499 vmax=180,
525 ax.ax_profile.grid(b=True, axis='x')
500 cmap=plt.get_cmap(self.colormap_phase)
526 ax.plt_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
501 )
527 color="k", linestyle="dashed", lw=2)[0]
528 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
529 else:
502 else:
530 ax.plt.set_array(z[n].T.ravel())
503 ax.plt.set_array(phase.T.ravel())
531 ax.plt_dop.set_data(mean[n], y)
504 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
532 if self.showprofile:
533 ax.plt_profile.set_data(self.data['rti'][self.max_time][n], y)
534 ax.plt_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
535
505
536 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
537 size=8)
538 self.saveTime = self.max_time
506 self.saveTime = self.max_time
539
507
540
508
509 class PlotSpectraMeanData(PlotSpectraData):
510 '''
511 Plot for Spectra and Mean
512 '''
513 CODE = 'spc_mean'
514 colormap = 'jro'
515
516
541 class PlotRTIData(PlotData):
517 class PlotRTIData(PlotData):
518 '''
519 Plot for RTI data
520 '''
542
521
543 CODE = 'rti'
522 CODE = 'rti'
544 colormap = 'jro'
523 colormap = 'jro'
545
524
546 def setup(self):
525 def setup(self):
526 self.xaxis = 'time'
547 self.ncols = 1
527 self.ncols = 1
548 self.nrows = self.dataOut.nChannels
528 self.nrows = len(self.data.channels)
549 self.width = 10
529 self.nplots = len(self.data.channels)
550 #TODO : arreglar la altura de la figura, esta hardcodeada.
551 #Se arreglo, testear!
552 if self.ind_plt_ch:
553 self.height = 3.2#*self.nrows if self.nrows<6 else 12
554 else:
555 self.height = 2.2*self.nrows if self.nrows<6 else 12
556
557 '''
558 if self.nrows==1:
559 self.height += 1
560 '''
561 self.ylabel = 'Range [Km]'
530 self.ylabel = 'Range [Km]'
562 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
531 self.cb_label = 'dB'
563
532 self.titles = ['{} Channel {}'.format(self.CODE.upper(), x) for x in range(self.nrows)]
564 '''
565 Logica:
566 1) Si la variable ind_plt_ch es True, va a crear mas de 1 figura
567 2) guardamos "Figures" en una lista y "axes" en otra, quizas se deberia guardar el
568 axis dentro de "Figures" como un diccionario.
569 '''
570 if self.ind_plt_ch is False: #standard mode
571
572 if self.figure is None: #solo para la priemra vez
573 self.figure = plt.figure(figsize=(self.width, self.height),
574 edgecolor='k',
575 facecolor='w')
576 else:
577 self.figure.clf()
578 self.axes = []
579
580
581 for n in range(self.nrows):
582 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
583 #ax = self.figure(n+1)
584 ax.firsttime = True
585 self.axes.append(ax)
586
587 else : #append one figure foreach channel in channelList
588 if self.figurelist == None:
589 self.figurelist = []
590 for n in range(self.nrows):
591 self.figure = plt.figure(figsize=(self.width, self.height),
592 edgecolor='k',
593 facecolor='w')
594 #add always one subplot
595 self.figurelist.append(self.figure)
596
597 else : # cada dia nuevo limpia el axes, pero mantiene el figure
598 for eachfigure in self.figurelist:
599 eachfigure.clf() # eliminaria todas las figuras de la lista?
600 self.axes = []
601
602 for eachfigure in self.figurelist:
603 ax = eachfigure.add_subplot(1,1,1) #solo 1 axis por figura
604 #ax = self.figure(n+1)
605 ax.firsttime = True
606 #Cada figura tiene un distinto puntero
607 self.axes.append(ax)
608 #plt.close(eachfigure)
609
610
533
611 def plot(self):
534 def plot(self):
535 self.x = self.data.times
536 self.y = self.data.heights
537 self.z = self.data[self.CODE]
538 self.z = numpy.ma.masked_invalid(self.z)
612
539
613 if self.ind_plt_ch is False: #standard mode
614 self.x = np.array(self.times)
615 self.y = self.dataOut.getHeiRange()
616 self.z = []
617
618 for ch in range(self.nrows):
619 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
620
621 self.z = np.array(self.z)
622 for n, ax in enumerate(self.axes):
540 for n, ax in enumerate(self.axes):
623 x, y, z = self.fill_gaps(*self.decimate())
541 x, y, z = self.fill_gaps(*self.decimate())
624 if self.xmin is None:
542 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
625 xmin = self.min_time
543 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
626 else:
627 xmin = fromtimestamp(int(self.xmin), self.min_time)
628 if self.xmax is None:
629 xmax = xmin + self.xrange*60*60
630 else:
631 xmax = xmin + (self.xmax - self.xmin) * 60 * 60
632 self.zmin = self.zmin if self.zmin else np.min(self.z)
633 self.zmax = self.zmax if self.zmax else np.max(self.z)
634 if ax.firsttime:
544 if ax.firsttime:
635 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
545 ax.plt = ax.pcolormesh(x, y, z[n].T,
636 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
637 plot = ax.pcolormesh(x, y, z[n].T,
638 vmin=self.zmin,
546 vmin=self.zmin,
639 vmax=self.zmax,
547 vmax=self.zmax,
640 cmap=plt.get_cmap(self.colormap)
548 cmap=plt.get_cmap(self.colormap)
641 )
549 )
642 divider = make_axes_locatable(ax)
550 if self.showprofile:
643 cax = divider.new_horizontal(size='2%', pad=0.05)
551 ax.plot_profile= self.pf_axes[n].plot(self.data['rti'][n][-1], self.y)[0]
644 self.figure.add_axes(cax)
552 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
645 plt.colorbar(plot, cax)
553 color="k", linestyle="dashed", lw=1)[0]
646 ax.set_ylim(self.ymin, self.ymax)
647 ax.xaxis.set_major_formatter(FuncFormatter(func))
648 ax.xaxis.set_major_locator(LinearLocator(6))
649 ax.set_ylabel(self.ylabel)
650 # if self.xmin is None:
651 # xmin = self.min_time
652 # else:
653 # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
654 # datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
655
656 ax.set_xlim(xmin, xmax)
657 ax.firsttime = False
658 else:
554 else:
659 ax.collections.remove(ax.collections[0])
555 ax.collections.remove(ax.collections[0])
660 ax.set_xlim(xmin, xmax)
556 ax.plt = ax.pcolormesh(x, y, z[n].T,
661 plot = ax.pcolormesh(x, y, z[n].T,
662 vmin=self.zmin,
663 vmax=self.zmax,
664 cmap=plt.get_cmap(self.colormap)
665 )
666 ax.set_title('{} {}'.format(self.titles[n],
667 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
668 size=8)
669
670 self.saveTime = self.min_time
671 else :
672 self.x = np.array(self.times)
673 self.y = self.dataOut.getHeiRange()
674 self.z = []
675
676 for ch in range(self.nrows):
677 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
678
679 self.z = np.array(self.z)
680 for n, eachfigure in enumerate(self.figurelist): #estaba ax in axes
681
682 x, y, z = self.fill_gaps(*self.decimate())
683 xmin = self.min_time
684 xmax = xmin+self.xrange*60*60
685 self.zmin = self.zmin if self.zmin else np.min(self.z)
686 self.zmax = self.zmax if self.zmax else np.max(self.z)
687 if self.axes[n].firsttime:
688 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
689 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
690 plot = self.axes[n].pcolormesh(x, y, z[n].T,
691 vmin=self.zmin,
692 vmax=self.zmax,
693 cmap=plt.get_cmap(self.colormap)
694 )
695 divider = make_axes_locatable(self.axes[n])
696 cax = divider.new_horizontal(size='2%', pad=0.05)
697 eachfigure.add_axes(cax)
698 #self.figure2.add_axes(cax)
699 plt.colorbar(plot, cax)
700 self.axes[n].set_ylim(self.ymin, self.ymax)
701
702 self.axes[n].xaxis.set_major_formatter(FuncFormatter(func))
703 self.axes[n].xaxis.set_major_locator(LinearLocator(6))
704
705 self.axes[n].set_ylabel(self.ylabel)
706
707 if self.xmin is None:
708 xmin = self.min_time
709 else:
710 xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
711 datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
712
713 self.axes[n].set_xlim(xmin, xmax)
714 self.axes[n].firsttime = False
715 else:
716 self.axes[n].collections.remove(self.axes[n].collections[0])
717 self.axes[n].set_xlim(xmin, xmax)
718 plot = self.axes[n].pcolormesh(x, y, z[n].T,
719 vmin=self.zmin,
557 vmin=self.zmin,
720 vmax=self.zmax,
558 vmax=self.zmax,
721 cmap=plt.get_cmap(self.colormap)
559 cmap=plt.get_cmap(self.colormap)
722 )
560 )
723 self.axes[n].set_title('{} {}'.format(self.titles[n],
561 if self.showprofile:
724 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
562 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
725 size=8)
563 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y)
726
564
727 self.saveTime = self.min_time
565 self.saveTime = self.min_time
728
566
729
567
730 class PlotCOHData(PlotRTIData):
568 class PlotCOHData(PlotRTIData):
569 '''
570 Plot for Coherence data
571 '''
731
572
732 CODE = 'coh'
573 CODE = 'coh'
733
574
734 def setup(self):
575 def setup(self):
735
576 self.xaxis = 'time'
736 self.ncols = 1
577 self.ncols = 1
737 self.nrows = self.dataOut.nPairs
578 self.nrows = len(self.data.pairs)
738 self.width = 10
579 self.nplots = len(self.data.pairs)
739 self.height = 2.2*self.nrows if self.nrows<6 else 12
740 self.ind_plt_ch = False #just for coherence and phase
741 if self.nrows==1:
742 self.height += 1
743 self.ylabel = 'Range [Km]'
580 self.ylabel = 'Range [Km]'
744 self.titles = ['{} Ch{} * Ch{}'.format(self.CODE.upper(), x[0], x[1]) for x in self.dataOut.pairsList]
581 if self.CODE == 'coh':
745
582 self.cb_label = ''
746 if self.figure is None:
583 self.titles = ['Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
747 self.figure = plt.figure(figsize=(self.width, self.height),
748 edgecolor='k',
749 facecolor='w')
750 else:
584 else:
751 self.figure.clf()
585 self.cb_label = 'Degrees'
752 self.axes = []
586 self.titles = ['Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
753
587
754 for n in range(self.nrows):
588
755 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
589 class PlotPHASEData(PlotCOHData):
756 ax.firsttime = True
590 '''
757 self.axes.append(ax)
591 Plot for Phase map data
592 '''
593
594 CODE = 'phase'
595 colormap = 'seismic'
758
596
759
597
760 class PlotNoiseData(PlotData):
598 class PlotNoiseData(PlotData):
599 '''
600 Plot for noise
601 '''
602
761 CODE = 'noise'
603 CODE = 'noise'
762
604
763 def setup(self):
605 def setup(self):
764
606 self.xaxis = 'time'
765 self.ncols = 1
607 self.ncols = 1
766 self.nrows = 1
608 self.nrows = 1
767 self.width = 10
609 self.nplots = 1
768 self.height = 3.2
769 self.ylabel = 'Intensity [dB]'
610 self.ylabel = 'Intensity [dB]'
770 self.titles = ['Noise']
611 self.titles = ['Noise']
771
612 self.colorbar = False
772 if self.figure is None:
773 self.figure = plt.figure(figsize=(self.width, self.height),
774 edgecolor='k',
775 facecolor='w')
776 else:
777 self.figure.clf()
778 self.axes = []
779
780 self.ax = self.figure.add_subplot(self.nrows, self.ncols, 1)
781 self.ax.firsttime = True
782
783 def plot(self):
784
785 x = self.times
786 xmin = self.min_time
787 xmax = xmin+self.xrange*60*60
788 if self.ax.firsttime:
789 for ch in self.dataOut.channelList:
790 y = [self.data[self.CODE][t][ch] for t in self.times]
791 self.ax.plot(x, y, lw=1, label='Ch{}'.format(ch))
792 self.ax.firsttime = False
793 self.ax.xaxis.set_major_formatter(FuncFormatter(func))
794 self.ax.xaxis.set_major_locator(LinearLocator(6))
795 self.ax.set_ylabel(self.ylabel)
796 plt.legend()
797 else:
798 for ch in self.dataOut.channelList:
799 y = [self.data[self.CODE][t][ch] for t in self.times]
800 self.ax.lines[ch].set_data(x, y)
801
802 self.ax.set_xlim(xmin, xmax)
803 self.ax.set_ylim(min(y)-5, max(y)+5)
804 self.saveTime = self.min_time
805
806
807 class PlotWindProfilerData(PlotRTIData):
808
809 CODE = 'wind'
810 colormap = 'seismic'
811
812 def setup(self):
813 self.ncols = 1
814 self.nrows = self.dataOut.data_output.shape[0]
815 self.width = 10
816 self.height = 2.2*self.nrows
817 self.ylabel = 'Height [Km]'
818 self.titles = ['Zonal Wind' ,'Meridional Wind', 'Vertical Wind']
819 self.clabels = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
820 self.windFactor = [1, 1, 100]
821
822 if self.figure is None:
823 self.figure = plt.figure(figsize=(self.width, self.height),
824 edgecolor='k',
825 facecolor='w')
826 else:
827 self.figure.clf()
828 self.axes = []
829
830 for n in range(self.nrows):
831 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
832 ax.firsttime = True
833 self.axes.append(ax)
834
613
835 def plot(self):
614 def plot(self):
836
615
837 self.x = np.array(self.times)
616 x = self.data.times
838 self.y = self.dataOut.heightList
839 self.z = []
840
841 for ch in range(self.nrows):
842 self.z.append([self.data['output'][t][ch] for t in self.times])
843
844 self.z = np.array(self.z)
845 self.z = numpy.ma.masked_invalid(self.z)
846
847 cmap=plt.get_cmap(self.colormap)
848 cmap.set_bad('black', 1.)
849
850 for n, ax in enumerate(self.axes):
851 x, y, z = self.fill_gaps(*self.decimate())
852 xmin = self.min_time
617 xmin = self.min_time
853 xmax = xmin+self.xrange*60*60
618 xmax = xmin+self.xrange*60*60
854 if ax.firsttime:
619 Y = self.data[self.CODE]
855 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
856 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
857 self.zmax = self.zmax if self.zmax else numpy.nanmax(abs(self.z[:-1, :]))
858 self.zmin = self.zmin if self.zmin else -self.zmax
859
860 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
861 vmin=self.zmin,
862 vmax=self.zmax,
863 cmap=cmap
864 )
865 divider = make_axes_locatable(ax)
866 cax = divider.new_horizontal(size='2%', pad=0.05)
867 self.figure.add_axes(cax)
868 cb = plt.colorbar(plot, cax)
869 cb.set_label(self.clabels[n])
870 ax.set_ylim(self.ymin, self.ymax)
871
620
872 ax.xaxis.set_major_formatter(FuncFormatter(func))
621 if self.axes[0].firsttime:
873 ax.xaxis.set_major_locator(LinearLocator(6))
622 for ch in self.data.channels:
874
623 y = Y[ch]
875 ax.set_ylabel(self.ylabel)
624 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
876
625 plt.legend()
877 ax.set_xlim(xmin, xmax)
878 ax.firsttime = False
879 else:
626 else:
880 ax.collections.remove(ax.collections[0])
627 for ch in self.data.channels:
881 ax.set_xlim(xmin, xmax)
628 y = Y[ch]
882 plot = ax.pcolormesh(x, y, z[n].T*self.windFactor[n],
629 self.axes[0].lines[ch].set_data(x, y)
883 vmin=self.zmin,
884 vmax=self.zmax,
885 cmap=plt.get_cmap(self.colormap)
886 )
887 ax.set_title('{} {}'.format(self.titles[n],
888 datetime.datetime.fromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
889 size=8)
890
630
631 self.ymin = numpy.nanmin(Y) - 5
632 self.ymax = numpy.nanmax(Y) + 5
891 self.saveTime = self.min_time
633 self.saveTime = self.min_time
892
634
893
635
894 class PlotSNRData(PlotRTIData):
636 class PlotSNRData(PlotRTIData):
637 '''
638 Plot for SNR Data
639 '''
640
895 CODE = 'snr'
641 CODE = 'snr'
896 colormap = 'jet'
642 colormap = 'jet'
897
643
644
898 class PlotDOPData(PlotRTIData):
645 class PlotDOPData(PlotRTIData):
646 '''
647 Plot for DOPPLER Data
648 '''
649
899 CODE = 'dop'
650 CODE = 'dop'
900 colormap = 'jet'
651 colormap = 'jet'
901
652
902
653
903 class PlotPHASEData(PlotCOHData):
904 CODE = 'phase'
905 colormap = 'seismic'
906
907
908 class PlotSkyMapData(PlotData):
654 class PlotSkyMapData(PlotData):
655 '''
656 Plot for meteors detection data
657 '''
909
658
910 CODE = 'met'
659 CODE = 'met'
911
660
@@ -932,7 +681,7 class PlotSkyMapData(PlotData):
932
681
933 def plot(self):
682 def plot(self):
934
683
935 arrayParameters = np.concatenate([self.data['param'][t] for t in self.times])
684 arrayParameters = numpy.concatenate([self.data['param'][t] for t in self.data.times])
936 error = arrayParameters[:,-1]
685 error = arrayParameters[:,-1]
937 indValid = numpy.where(error == 0)[0]
686 indValid = numpy.where(error == 0)[0]
938 finalMeteor = arrayParameters[indValid,:]
687 finalMeteor = arrayParameters[indValid,:]
@@ -962,3 +711,72 class PlotSkyMapData(PlotData):
962 self.ax.set_title(title, size=8)
711 self.ax.set_title(title, size=8)
963
712
964 self.saveTime = self.max_time
713 self.saveTime = self.max_time
714
715 class PlotParamData(PlotRTIData):
716 '''
717 Plot for data_param object
718 '''
719
720 CODE = 'param'
721 colormap = 'seismic'
722
723 def setup(self):
724 self.xaxis = 'time'
725 self.ncols = 1
726 self.nrows = self.data.shape(self.CODE)[0]
727 self.nplots = self.nrows
728 if self.showSNR:
729 self.nrows += 1
730
731 self.ylabel = 'Height [Km]'
732 self.titles = self.data.parameters \
733 if self.data.parameters else ['Param {}'.format(x) for x in xrange(self.nrows)]
734 if self.showSNR:
735 self.titles.append('SNR')
736
737 def plot(self):
738 self.data.normalize_heights()
739 self.x = self.data.times
740 self.y = self.data.heights
741 if self.showSNR:
742 self.z = numpy.concatenate(
743 (self.data[self.CODE], self.data['snr'])
744 )
745 else:
746 self.z = self.data[self.CODE]
747
748 self.z = numpy.ma.masked_invalid(self.z)
749
750 for n, ax in enumerate(self.axes):
751
752 x, y, z = self.fill_gaps(*self.decimate())
753
754 if ax.firsttime:
755 if self.zlimits is not None:
756 self.zmin, self.zmax = self.zlimits[n]
757 self.zmax = self.zmax if self.zmax is not None else numpy.nanmax(abs(self.z[:-1, :]))
758 self.zmin = self.zmin if self.zmin is not None else -self.zmax
759 ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n],
760 vmin=self.zmin,
761 vmax=self.zmax,
762 cmap=self.cmaps[n]
763 )
764 else:
765 if self.zlimits is not None:
766 self.zmin, self.zmax = self.zlimits[n]
767 ax.collections.remove(ax.collections[0])
768 ax.plt = ax.pcolormesh(x, y, z[n, :, :].T*self.factors[n],
769 vmin=self.zmin,
770 vmax=self.zmax,
771 cmap=self.cmaps[n]
772 )
773
774 self.saveTime = self.min_time
775
776 class PlotOuputData(PlotParamData):
777 '''
778 Plot data_output object
779 '''
780
781 CODE = 'output'
782 colormap = 'seismic' No newline at end of file
Show file before | Show file after | Hide comments
@@ -267,7 +267,7 class AMISRReader(ProcessingUnit):
267 self.dirnameList = new_dirnameList
267 self.dirnameList = new_dirnameList
268 return 1
268 return 1
269
269
270 def __searchFilesOnline(self,
270 def searchFilesOnLine(self,
271 path,
271 path,
272 walk=True):
272 walk=True):
273
273
@@ -287,7 +287,7 class AMISRReader(ProcessingUnit):
287 return
287 return
288
288
289
289
290 def __searchFilesOffline(self,
290 def searchFilesOffLine(self,
291 path,
291 path,
292 startDate,
292 startDate,
293 endDate,
293 endDate,
@@ -494,9 +494,9 class AMISRReader(ProcessingUnit):
494 self.online = online
494 self.online = online
495 if not(online):
495 if not(online):
496 #Busqueda de archivos offline
496 #Busqueda de archivos offline
497 self.__searchFilesOffline(path, startDate, endDate, startTime, endTime, walk)
497 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
498 else:
498 else:
499 self.__searchFilesOnline(path, walk)
499 self.searchFilesOnLine(path, walk)
500
500
501 if not(self.filenameList):
501 if not(self.filenameList):
502 print "There is no files into the folder: %s"%(path)
502 print "There is no files into the folder: %s"%(path)
Show file before | Show file after | Hide comments
@@ -542,7 +542,6 class JRODataIO:
542
542
543 class JRODataReader(JRODataIO):
543 class JRODataReader(JRODataIO):
544
544
545 firstTime = True
546 online = 0
545 online = 0
547
546
548 realtime = 0
547 realtime = 0
@@ -579,7 +578,6 class JRODataReader(JRODataIO):
579
578
580 selBlocktime = None
579 selBlocktime = None
581
580
582 onlineWithDate = False
583 def __init__(self):
581 def __init__(self):
584
582
585 """
583 """
@@ -603,7 +601,7 class JRODataReader(JRODataIO):
603
601
604 raise NotImplementedError
602 raise NotImplementedError
605
603
606 def __searchFilesOffLine(self,
604 def searchFilesOffLine(self,
607 path,
605 path,
608 startDate=None,
606 startDate=None,
609 endDate=None,
607 endDate=None,
@@ -612,10 +610,10 class JRODataReader(JRODataIO):
612 set=None,
610 set=None,
613 expLabel='',
611 expLabel='',
614 ext='.r',
612 ext='.r',
615 queue=None,
616 cursor=None,
613 cursor=None,
617 skip=None,
614 skip=None,
618 walk=True):
615 walk=True):
616
619 self.filenameList = []
617 self.filenameList = []
620 self.datetimeList = []
618 self.datetimeList = []
621
619
@@ -624,8 +622,7 class JRODataReader(JRODataIO):
624 dateList, pathList = self.findDatafiles(path, startDate, endDate, expLabel, ext, walk, include_path=True)
622 dateList, pathList = self.findDatafiles(path, startDate, endDate, expLabel, ext, walk, include_path=True)
625
623
626 if dateList == []:
624 if dateList == []:
627 # print "[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" %(startDate, endDate, ext, path)
625 return [], []
628 return None, None
629
626
630 if len(dateList) > 1:
627 if len(dateList) > 1:
631 print "[Reading] Data found for date range [%s - %s]: total days = %d" %(startDate, endDate, len(dateList))
628 print "[Reading] Data found for date range [%s - %s]: total days = %d" %(startDate, endDate, len(dateList))
@@ -636,18 +633,15 class JRODataReader(JRODataIO):
636 datetimeList = []
633 datetimeList = []
637
634
638 for thisPath in pathList:
635 for thisPath in pathList:
639 # thisPath = pathList[pathDict[file]]
640
636
641 fileList = glob.glob1(thisPath, "*%s" %ext)
637 fileList = glob.glob1(thisPath, "*%s" %ext)
642 fileList.sort()
638 fileList.sort()
643
639
644 skippedFileList = []
640 skippedFileList = []
645
641
646 if cursor is not None and skip is not None:
642 if cursor is not None andk skip is not None:
647 # if cursor*skip > len(fileList):
643
648 if skip == 0:
644 if skip == 0:
649 if queue is not None:
650 queue.put(len(fileList))
651 skippedFileList = []
645 skippedFileList = []
652 else:
646 else:
653 skippedFileList = fileList[cursor*skip: cursor*skip + skip]
647 skippedFileList = fileList[cursor*skip: cursor*skip + skip]
@@ -672,19 +666,20 class JRODataReader(JRODataIO):
672
666
673 if not(filenameList):
667 if not(filenameList):
674 print "[Reading] Time range selected invalid [%s - %s]: No *%s files in %s)" %(startTime, endTime, ext, path)
668 print "[Reading] Time range selected invalid [%s - %s]: No *%s files in %s)" %(startTime, endTime, ext, path)
675 return None, None
669 return [], []
676
670
677 print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)
671 print "[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)
678 print
672 print
679
673
680 for i in range(len(filenameList)):
674 # for i in range(len(filenameList)):
681 print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
675 # print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
682
676
683 self.filenameList = filenameList
677 self.filenameList = filenameList
684 self.datetimeList = datetimeList
678 self.datetimeList = datetimeList
679
685 return pathList, filenameList
680 return pathList, filenameList
686
681
687 def __searchFilesOnLine(self, path, expLabel="", ext=None, walk=True, set=None, startDate=None, startTime=None):
682 def __searchFilesOnLine(self, path, expLabel = "", ext = None, walk=True, set=None):
688
683
689 """
684 """
690 Busca el ultimo archivo de la ultima carpeta (determinada o no por startDateTime) y
685 Busca el ultimo archivo de la ultima carpeta (determinada o no por startDateTime) y
@@ -697,7 +692,7 class JRODataReader(JRODataIO):
697
692
698 ext : extension de los files
693 ext : extension de los files
699
694
700 walk : Si es habilitado no realiza busquedas dentro de los subdirectorios (doypath)
695 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
701
696
702 Return:
697 Return:
703 directory : eL directorio donde esta el file encontrado
698 directory : eL directorio donde esta el file encontrado
@@ -708,8 +703,6 class JRODataReader(JRODataIO):
708
703
709
704
710 """
705 """
711 pathList = None
712 filenameList = None
713 if not os.path.isdir(path):
706 if not os.path.isdir(path):
714 return None, None, None, None, None, None
707 return None, None, None, None, None, None
715
708
@@ -810,7 +803,6 class JRODataReader(JRODataIO):
810 Excepciones:
803 Excepciones:
811 Si un determinado file no puede ser abierto
804 Si un determinado file no puede ser abierto
812 """
805 """
813
814 nFiles = 0
806 nFiles = 0
815 fileOk_flag = False
807 fileOk_flag = False
816 firstTime_flag = True
808 firstTime_flag = True
@@ -883,34 +875,13 class JRODataReader(JRODataIO):
883 def setNextFile(self):
875 def setNextFile(self):
884 if self.fp != None:
876 if self.fp != None:
885 self.fp.close()
877 self.fp.close()
878
886 if self.online:
879 if self.online:
887 newFile = self.__setNextFileOnline()
880 newFile = self.__setNextFileOnline()
888 else:
881 else:
889 newFile = self.__setNextFileOffline()
882 newFile = self.__setNextFileOffline()
883
890 if not(newFile):
884 if not(newFile):
891 if self.onlineWithDate is True:
892 self.onlineWithDate=False
893 self.online = True
894 self.firstTime = False
895 self.setup(
896 path=self.path,
897 startDate=self.startDate,
898 endDate=self.endDate,
899 startTime=self.startTime ,
900 endTime=self.endTime,
901 set=self.set,
902 expLabel=self.expLabel,
903 ext=self.ext,
904 online=self.online,
905 delay=self.delay,
906 walk=self.walk,
907 getblock=self.getblock,
908 nTxs=self.nTxs,
909 realtime=self.realtime,
910 blocksize=self.blocksize,
911 blocktime=self.blocktime
912 )
913 return 1
914 print '[Reading] No more files to read'
885 print '[Reading] No more files to read'
915 return 0
886 return 0
916
887
@@ -1066,11 +1037,13 class JRODataReader(JRODataIO):
1066 #Skip block out of startTime and endTime
1037 #Skip block out of startTime and endTime
1067 while True:
1038 while True:
1068 if not(self.__setNewBlock()):
1039 if not(self.__setNewBlock()):
1069 print 'returning'
1070 return 0
1040 return 0
1041
1071 if not(self.readBlock()):
1042 if not(self.readBlock()):
1072 return 0
1043 return 0
1044
1073 self.getBasicHeader()
1045 self.getBasicHeader()
1046
1074 if not isTimeInRange(self.dataOut.datatime.time(), self.startTime, self.endTime):
1047 if not isTimeInRange(self.dataOut.datatime.time(), self.startTime, self.endTime):
1075
1048
1076 print "[Reading] Block No. %d/%d -> %s [Skipping]" %(self.nReadBlocks,
1049 print "[Reading] Block No. %d/%d -> %s [Skipping]" %(self.nReadBlocks,
@@ -1292,66 +1265,40 class JRODataReader(JRODataIO):
1292 realtime=False,
1265 realtime=False,
1293 blocksize=None,
1266 blocksize=None,
1294 blocktime=None,
1267 blocktime=None,
1268 skip=None,
1269 cursor=None,
1270 warnings=True,
1295 verbose=True,
1271 verbose=True,
1272 server=None,
1296 **kwargs):
1273 **kwargs):
1297
1274 if server is not None:
1275 if 'tcp://' in server:
1276 address = server
1277 else:
1278 address = 'ipc:///tmp/%s' % server
1279 self.server = address
1280 self.context = zmq.Context()
1281 self.receiver = self.context.socket(zmq.PULL)
1282 self.receiver.connect(self.server)
1283 time.sleep(0.5)
1284 print '[Starting] ReceiverData from {}'.format(self.server)
1285 else:
1286 self.server = None
1298 if path == None:
1287 if path == None:
1299 raise ValueError, "[Reading] The path is not valid"
1288 raise ValueError, "[Reading] The path is not valid"
1300
1289
1301
1302 if ext == None:
1290 if ext == None:
1303 ext = self.ext
1291 ext = self.ext
1304
1292
1305 self.verbose=verbose
1306 self.path = path
1307 self.startDate = startDate
1308 self.endDate = endDate
1309 self.startTime = startTime
1310 self.endTime = endTime
1311 self.set = set
1312 self.expLabel = expLabel
1313 self.ext = ext
1314 self.online = online
1315 self.delay = delay
1316 self.walk = walk
1317 self.getblock = getblock
1318 self.nTxs = nTxs
1319 self.realtime = realtime
1320 self.blocksize = blocksize
1321 self.blocktime = blocktime
1322
1323
1324 if self.firstTime is True:
1325 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
1326 startTime=startTime, endTime=endTime,
1327 set=set, expLabel=expLabel, ext=ext,
1328 walk=walk)
1329 if filenameList is not None: filenameList = filenameList[:-1]
1330
1331 if pathList is not None and filenameList is not None and online:
1332 self.onlineWithDate = True
1333 online = False
1334 self.fileIndex = -1
1335 self.pathList = pathList
1336 self.filenameList = filenameList
1337 file_name = os.path.basename(filenameList[-1])
1338 basename, ext = os.path.splitext(file_name)
1339 last_set = int(basename[-3:])
1340
1341 if online:
1293 if online:
1342 print "[Reading] Searching files in online mode..."
1294 print "[Reading] Searching files in online mode..."
1343
1295
1344 for nTries in range(self.nTries):
1296 for nTries in range( self.nTries ):
1345 fullpath, foldercounter, file, year, doy, set = self.__searchFilesOnLine(path=path,
1297 fullpath, foldercounter, file, year, doy, set = self.__searchFilesOnLine(path=path, expLabel=expLabel, ext=ext, walk=walk, set=set)
1346 expLabel=expLabel,
1347 ext=ext,
1348 walk=walk,
1349 startDate=startDate,
1350 startTime=startTime,
1351 set=set)
1352
1298
1353 if fullpath:
1299 if fullpath:
1354 break
1300 break
1301
1355 print '[Reading] Waiting %0.2f sec for an valid file in %s: try %02d ...' % (self.delay, path, nTries+1)
1302 print '[Reading] Waiting %0.2f sec for an valid file in %s: try %02d ...' % (self.delay, path, nTries+1)
1356 sleep( self.delay )
1303 sleep( self.delay )
1357
1304
@@ -1367,18 +1314,13 class JRODataReader(JRODataIO):
1367 last_set = None
1314 last_set = None
1368 else:
1315 else:
1369 print "[Reading] Searching files in offline mode ..."
1316 print "[Reading] Searching files in offline mode ..."
1370 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
1317 pathList, filenameList = self.searchFilesOffLine(path, startDate=startDate, endDate=endDate,
1371 startTime=startTime, endTime=endTime,
1318 startTime=startTime, endTime=endTime,
1372 set=set, expLabel=expLabel, ext=ext,
1319 set=set, expLabel=expLabel, ext=ext,
1373 walk=walk)
1320 walk=walk, cursor=cursor,
1321 skip=skip)
1374
1322
1375 if not(pathList):
1323 if not(pathList):
1376 # print "[Reading] No *%s files in %s (%s - %s)"%(ext, path,
1377 # datetime.datetime.combine(startDate,startTime).ctime(),
1378 # datetime.datetime.combine(endDate,endTime).ctime())
1379
1380 # sys.exit(-1)
1381
1382 self.fileIndex = -1
1324 self.fileIndex = -1
1383 self.pathList = []
1325 self.pathList = []
1384 self.filenameList = []
1326 self.filenameList = []
@@ -1391,7 +1333,6 class JRODataReader(JRODataIO):
1391 basename, ext = os.path.splitext(file_name)
1333 basename, ext = os.path.splitext(file_name)
1392 last_set = int(basename[-3:])
1334 last_set = int(basename[-3:])
1393
1335
1394
1395 self.online = online
1336 self.online = online
1396 self.realtime = realtime
1337 self.realtime = realtime
1397 self.delay = delay
1338 self.delay = delay
@@ -1402,11 +1343,13 class JRODataReader(JRODataIO):
1402 self.startTime = startTime
1343 self.startTime = startTime
1403 self.endTime = endTime
1344 self.endTime = endTime
1404
1345
1405
1406 #Added-----------------
1346 #Added-----------------
1407 self.selBlocksize = blocksize
1347 self.selBlocksize = blocksize
1408 self.selBlocktime = blocktime
1348 self.selBlocktime = blocktime
1409
1349
1350 # Verbose-----------
1351 self.verbose = verbose
1352 self.warnings = warnings
1410
1353
1411 if not(self.setNextFile()):
1354 if not(self.setNextFile()):
1412 if (startDate!=None) and (endDate!=None):
1355 if (startDate!=None) and (endDate!=None):
@@ -1517,7 +1460,6 class JRODataReader(JRODataIO):
1517 server=None,
1460 server=None,
1518 verbose=True, **kwargs):
1461 verbose=True, **kwargs):
1519 if not(self.isConfig):
1462 if not(self.isConfig):
1520 # self.dataOut = dataOut
1521 self.setup( path=path,
1463 self.setup(path=path,
1522 startDate=startDate,
1464 startDate=startDate,
1523 endDate=endDate,
1465 endDate=endDate,
@@ -1534,12 +1476,11 class JRODataReader(JRODataIO):
1534 realtime=realtime,
1476 realtime=realtime,
1535 blocksize=blocksize,
1477 blocksize=blocksize,
1536 blocktime=blocktime,
1478 blocktime=blocktime,
1537 queue=queue,
1538 skip=skip,
1479 skip=skip,
1539 cursor=cursor,
1480 cursor=cursor,
1540 warnings=warnings,
1481 warnings=warnings,
1541 server=server,
1482 server=server,
1542 verbose=verbose, **kwargs)
1483 verbose=verbose)
1543 self.isConfig = True
1484 self.isConfig = True
1544 if server is None:
1485 if server is None:
1545 self.getData()
1486 self.getData()
@@ -1791,7 +1732,7 class JRODataWriter(JRODataIO):
1791
1732
1792 return 1
1733 return 1
1793
1734
1794 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, verbose=True):
1735 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1795 """
1736 """
1796 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1737 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1797
1738
Show file before | Show file after | Hide comments
@@ -453,7 +453,7 class FitsReader(ProcessingUnit):
453 # self.blockIndex = 1
453 # self.blockIndex = 1
454 return 1
454 return 1
455
455
456 def __searchFilesOffLine(self,
456 def searchFilesOffLine(self,
457 path,
457 path,
458 startDate,
458 startDate,
459 endDate,
459 endDate,
@@ -559,7 +559,7 class FitsReader(ProcessingUnit):
559
559
560 if not(online):
560 if not(online):
561 print "Searching files in offline mode ..."
561 print "Searching files in offline mode ..."
562 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
562 pathList, filenameList = self.searchFilesOffLine(path, startDate=startDate, endDate=endDate,
563 startTime=startTime, endTime=endTime,
563 startTime=startTime, endTime=endTime,
564 set=set, expLabel=expLabel, ext=ext,
564 set=set, expLabel=expLabel, ext=ext,
565 walk=walk)
565 walk=walk)
Show file before | Show file after | Hide comments
@@ -415,7 +415,7 class HFReader(ProcessingUnit):
415
415
416
416
417
417
418 def __searchFilesOffline(self,
418 def searchFilesOffLine(self,
419 path,
419 path,
420 startDate,
420 startDate,
421 endDate,
421 endDate,
@@ -438,7 +438,7 class HFReader(ProcessingUnit):
438
438
439 return
439 return
440
440
441 def __searchFilesOnline(self,
441 def searchFilesOnLine(self,
442 path,
442 path,
443 expLabel= "",
443 expLabel= "",
444 ext=None,
444 ext=None,
@@ -636,10 +636,10 class HFReader(ProcessingUnit):
636 if not(online):
636 if not(online):
637 print "Searching files in offline mode..."
637 print "Searching files in offline mode..."
638
638
639 self.__searchFilesOffline(path, startDate, endDate, ext, startTime, endTime, walk)
639 self.searchFilesOffLine(path, startDate, endDate, ext, startTime, endTime, walk)
640 else:
640 else:
641 print "Searching files in online mode..."
641 print "Searching files in online mode..."
642 self.__searchFilesOnline(path, walk,ext,set=set)
642 self.searchFilesOnLine(path, walk,ext,set=set)
643 if set==None:
643 if set==None:
644 pass
644 pass
645 else:
645 else:
@@ -647,7 +647,7 class HFReader(ProcessingUnit):
647
647
648 # for nTries in range(self.nTries):
648 # for nTries in range(self.nTries):
649 #
649 #
650 # fullpath,file,year,month,day,set = self.__searchFilesOnline(path=path,expLabel=expLabel,ext=ext, walk=walk,set=set)
650 # fullpath,file,year,month,day,set = self.searchFilesOnLine(path=path,expLabel=expLabel,ext=ext, walk=walk,set=set)
651 #
651 #
652 # if fullpath:
652 # if fullpath:
653 # break
653 # break
Show file before | Show file after | Hide comments
@@ -106,9 +106,9 class AMISRReader(ProcessingUnit):
106 #self.findFiles()
106 #self.findFiles()
107 if not(online):
107 if not(online):
108 #Busqueda de archivos offline
108 #Busqueda de archivos offline
109 self.__searchFilesOffline(path, startDate, endDate, startTime, endTime, walk)
109 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
110 else:
110 else:
111 self.__searchFilesOnline(path, startDate, endDate, startTime,endTime,walk)
111 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
112
112
113 if not(self.filenameList):
113 if not(self.filenameList):
114 print "There is no files into the folder: %s"%(path)
114 print "There is no files into the folder: %s"%(path)
@@ -329,7 +329,7 class AMISRReader(ProcessingUnit):
329 self.dirnameList = new_dirnameList
329 self.dirnameList = new_dirnameList
330 return 1
330 return 1
331
331
332 def __searchFilesOnline(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
332 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
333 endTime=datetime.time(23,59,59),walk=True):
333 endTime=datetime.time(23,59,59),walk=True):
334
334
335 if endDate ==None:
335 if endDate ==None:
@@ -349,7 +349,7 class AMISRReader(ProcessingUnit):
349 return
349 return
350
350
351
351
352 def __searchFilesOffline(self,
352 def searchFilesOffLine(self,
353 path,
353 path,
354 startDate,
354 startDate,
355 endDate,
355 endDate,
Show file before | Show file after | Hide comments
@@ -97,7 +97,7 class ParamReader(ProcessingUnit):
97 self.timezone = 'lt'
97 self.timezone = 'lt'
98
98
99 print "[Reading] Searching files in offline mode ..."
99 print "[Reading] Searching files in offline mode ..."
100 pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate,
100 pathList, filenameList = self.searchFilesOffLine(path, startDate=startDate, endDate=endDate,
101 startTime=startTime, endTime=endTime,
101 startTime=startTime, endTime=endTime,
102 ext=ext, walk=walk)
102 ext=ext, walk=walk)
103
103
@@ -115,7 +115,7 class ParamReader(ProcessingUnit):
115
115
116 return
116 return
117
117
118 def __searchFilesOffLine(self,
118 def searchFilesOffLine(self,
119 path,
119 path,
120 startDate=None,
120 startDate=None,
121 endDate=None,
121 endDate=None,
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1338 lines changed) Show them Hide them
@@ -1,16 +1,56
1 import numpy
1 import numpy
2 import math
2 import math
3 from scipy import optimize, interpolate, signal, stats, ndimage
3 from scipy import optimize, interpolate, signal, stats, ndimage
4 import scipy
4 import re
5 import re
5 import datetime
6 import datetime
6 import copy
7 import copy
7 import sys
8 import sys
8 import importlib
9 import importlib
9 import itertools
10 import itertools
10
11 from multiprocessing import Pool, TimeoutError
12 from multiprocessing.pool import ThreadPool
13 import copy_reg
14 import cPickle
15 import types
16 from functools import partial
17 import time
18 #from sklearn.cluster import KMeans
19
20 import matplotlib.pyplot as plt
21
22 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
11 from jroproc_base import ProcessingUnit, Operation
23 from jroproc_base import ProcessingUnit, Operation
12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
24 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
25 from scipy import asarray as ar,exp
26 from scipy.optimize import curve_fit
27
28 import warnings
29 from numpy import NaN
30 from scipy.optimize.optimize import OptimizeWarning
31 warnings.filterwarnings('ignore')
32
33
34 SPEED_OF_LIGHT = 299792458
35
36
37 '''solving pickling issue'''
13
38
39 def _pickle_method(method):
40 func_name = method.im_func.__name__
41 obj = method.im_self
42 cls = method.im_class
43 return _unpickle_method, (func_name, obj, cls)
44
45 def _unpickle_method(func_name, obj, cls):
46 for cls in cls.mro():
47 try:
48 func = cls.__dict__[func_name]
49 except KeyError:
50 pass
51 else:
52 break
53 return func.__get__(obj, cls)
14
54
15 class ParametersProc(ProcessingUnit):
55 class ParametersProc(ProcessingUnit):
16
56
@@ -83,12 +123,32 class ParametersProc(ProcessingUnit):
83 self.dataOut.nIncohInt = self.dataIn.nIncohInt
123 self.dataOut.nIncohInt = self.dataIn.nIncohInt
84 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
124 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
85 self.dataOut.ippFactor = self.dataIn.ippFactor
125 self.dataOut.ippFactor = self.dataIn.ippFactor
86 #self.dataOut.normFactor = self.dataIn.getNormFactor()
126 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
127 self.dataOut.spc_noise = self.dataIn.getNoise()
128 self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
87 self.dataOut.pairsList = self.dataIn.pairsList
129 self.dataOut.pairsList = self.dataIn.pairsList
88 self.dataOut.groupList = self.dataIn.pairsList
130 self.dataOut.groupList = self.dataIn.pairsList
89 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
90 self.dataOut.flagNoData = False
131 self.dataOut.flagNoData = False
91
132
133 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
134 self.dataOut.ChanDist = self.dataIn.ChanDist
135 else: self.dataOut.ChanDist = None
136
137 if hasattr(self.dataIn, 'VelRange'): #Velocities range
138 self.dataOut.VelRange = self.dataIn.VelRange
139 else: self.dataOut.VelRange = None
140
141 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
142 self.dataOut.RadarConst = self.dataIn.RadarConst
143
144 if hasattr(self.dataIn, 'NPW'): #NPW
145 self.dataOut.NPW = self.dataIn.NPW
146
147 if hasattr(self.dataIn, 'COFA'): #COFA
148 self.dataOut.COFA = self.dataIn.COFA
149
150
151
92 #---------------------- Correlation Data ---------------------------
152 #---------------------- Correlation Data ---------------------------
93
153
94 if self.dataIn.type == "Correlation":
154 if self.dataIn.type == "Correlation":
@@ -108,7 +168,6 class ParametersProc(ProcessingUnit):
108
168
109 if self.dataIn.type == "Parameters":
169 if self.dataIn.type == "Parameters":
110 self.dataOut.copy(self.dataIn)
170 self.dataOut.copy(self.dataIn)
111 self.dataOut.utctimeInit = self.dataIn.utctime
112 self.dataOut.flagNoData = False
171 self.dataOut.flagNoData = False
113
172
114 return True
173 return True
@@ -119,6 +178,1186 class ParametersProc(ProcessingUnit):
119
178
120 return
179 return
121
180
181
182 def target(tups):
183
184 obj, args = tups
185 #print 'TARGETTT', obj, args
186 return obj.FitGau(args)
187
188 class GaussianFit(Operation):
189
190 '''
191 Function that fit of one and two generalized gaussians (gg) based
192 on the PSD shape across an "power band" identified from a cumsum of
193 the measured spectrum - noise.
194
195 Input:
196 self.dataOut.data_pre : SelfSpectra
197
198 Output:
199 self.dataOut.GauSPC : SPC_ch1, SPC_ch2
200
201 '''
202 def __init__(self, **kwargs):
203 Operation.__init__(self, **kwargs)
204 self.i=0
205
206
207 def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
208 """This routine will find a couple of generalized Gaussians to a power spectrum
209 input: spc
210 output:
211 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
212 """
213
214 self.spc = dataOut.data_pre[0].copy()
215
216
217 print 'SelfSpectra Shape', numpy.asarray(self.spc).shape
218
219
220 #plt.figure(50)
221 #plt.subplot(121)
222 #plt.plot(self.spc,'k',label='spc(66)')
223 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
224 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
225 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
226 #plt.legend()
227 #plt.title('DATOS A ALTURA DE 7500 METROS')
228 #plt.show()
229
230 self.Num_Hei = self.spc.shape[2]
231 #self.Num_Bin = len(self.spc)
232 self.Num_Bin = self.spc.shape[1]
233 self.Num_Chn = self.spc.shape[0]
234
235 Vrange = dataOut.abscissaList
236
237 #print 'self.spc2', numpy.asarray(self.spc).shape
238
239 GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
240 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
241 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
242 SPC_ch1[:] = numpy.NaN
243 SPC_ch2[:] = numpy.NaN
244
245
246 start_time = time.time()
247
248 noise_ = dataOut.spc_noise[0].copy()
249
250
251
252 pool = Pool(processes=self.Num_Chn)
253 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
254 objs = [self for __ in range(self.Num_Chn)]
255 attrs = zip(objs, args)
256 gauSPC = pool.map(target, attrs)
257 dataOut.GauSPC = numpy.asarray(gauSPC)
258 # ret = []
259 # for n in range(self.Num_Chn):
260 # self.FitGau(args[n])
261 # dataOut.GauSPC = ret
262
263
264
265 # for ch in range(self.Num_Chn):
266 #
267 # for ht in range(self.Num_Hei):
268 # #print (numpy.asarray(self.spc).shape)
269 # spc = numpy.asarray(self.spc)[ch,:,ht]
270 #
271 # #############################################
272 # # normalizing spc and noise
273 # # This part differs from gg1
274 # spc_norm_max = max(spc)
275 # spc = spc / spc_norm_max
276 # pnoise = pnoise / spc_norm_max
277 # #############################################
278 #
279 # if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
280 # fatspectra=1.0
281 # else:
282 # fatspectra=0.5
283 #
284 # wnoise = noise_ / spc_norm_max
285 # #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
286 # #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
287 # #if wnoise>1.1*pnoise: # to be tested later
288 # # wnoise=pnoise
289 # noisebl=wnoise*0.9; noisebh=wnoise*1.1
290 # spc=spc-wnoise
291 #
292 # minx=numpy.argmin(spc)
293 # spcs=numpy.roll(spc,-minx)
294 # cum=numpy.cumsum(spcs)
295 # tot_noise=wnoise * self.Num_Bin #64;
296 # #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
297 # #snr=tot_signal/tot_noise
298 # #snr=cum[-1]/tot_noise
299 #
300 # #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
301 #
302 # snr = sum(spcs)/tot_noise
303 # snrdB=10.*numpy.log10(snr)
304 #
305 # #if snrdB < -9 :
306 # # snrdB = numpy.NaN
307 # # continue
308 #
309 # #print 'snr',snrdB # , sum(spcs) , tot_noise
310 #
311 #
312 # #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
313 # # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
314 #
315 # cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
316 # cumlo=cummax*epsi;
317 # cumhi=cummax*(1-epsi)
318 # powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
319 #
320 # #if len(powerindex)==1:
321 # ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
322 # #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
323 # #elif len(powerindex)<4*fatspectra:
324 # #return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
325 #
326 # if len(powerindex) < 1:# case for powerindex 0
327 # continue
328 # powerlo=powerindex[0]
329 # powerhi=powerindex[-1]
330 # powerwidth=powerhi-powerlo
331 #
332 # firstpeak=powerlo+powerwidth/10.# first gaussian energy location
333 # secondpeak=powerhi-powerwidth/10.#second gaussian energy location
334 # midpeak=(firstpeak+secondpeak)/2.
335 # firstamp=spcs[int(firstpeak)]
336 # secondamp=spcs[int(secondpeak)]
337 # midamp=spcs[int(midpeak)]
338 # #x=numpy.spc.shape[1]
339 #
340 # #x=numpy.arange(64)
341 # x=numpy.arange( self.Num_Bin )
342 # y_data=spc+wnoise
343 #
344 # # single gaussian
345 # #shift0=numpy.mod(midpeak+minx,64)
346 # shift0=numpy.mod(midpeak+minx, self.Num_Bin )
347 # width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
348 # power0=2.
349 # amplitude0=midamp
350 # state0=[shift0,width0,amplitude0,power0,wnoise]
351 # #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
352 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
353 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
354 # # bnds = range of fft, power width, amplitude, power, noise
355 # lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
356 #
357 # chiSq1=lsq1[1];
358 # jack1= self.y_jacobian1(x,lsq1[0])
359 #
360 #
361 # try:
362 # sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
363 # except:
364 # std1=32.; sigmas1=numpy.ones(5)
365 # else:
366 # std1=sigmas1[0]
367 #
368 #
369 # if fatspectra<1.0 and powerwidth<4:
370 # choice=0
371 # Amplitude0=lsq1[0][2]
372 # shift0=lsq1[0][0]
373 # width0=lsq1[0][1]
374 # p0=lsq1[0][3]
375 # Amplitude1=0.
376 # shift1=0.
377 # width1=0.
378 # p1=0.
379 # noise=lsq1[0][4]
380 # #return (numpy.array([shift0,width0,Amplitude0,p0]),
381 # # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
382 #
383 # # two gaussians
384 # #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
385 # shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
386 # shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
387 # width0=powerwidth/6.;
388 # width1=width0
389 # power0=2.;
390 # power1=power0
391 # amplitude0=firstamp;
392 # amplitude1=secondamp
393 # state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
394 # #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
395 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
396 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
397 #
398 # lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
399 #
400 #
401 # chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
402 #
403 #
404 # try:
405 # sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
406 # except:
407 # std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
408 # else:
409 # std2a=sigmas2[0]; std2b=sigmas2[4]
410 #
411 #
412 #
413 # oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
414 #
415 # if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
416 # if oneG:
417 # choice=0
418 # else:
419 # w1=lsq2[0][1]; w2=lsq2[0][5]
420 # a1=lsq2[0][2]; a2=lsq2[0][6]
421 # p1=lsq2[0][3]; p2=lsq2[0][7]
422 # s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
423 # gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
424 #
425 # if gp1>gp2:
426 # if a1>0.7*a2:
427 # choice=1
428 # else:
429 # choice=2
430 # elif gp2>gp1:
431 # if a2>0.7*a1:
432 # choice=2
433 # else:
434 # choice=1
435 # else:
436 # choice=numpy.argmax([a1,a2])+1
437 # #else:
438 # #choice=argmin([std2a,std2b])+1
439 #
440 # else: # with low SNR go to the most energetic peak
441 # choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
442 #
443 # #print 'choice',choice
444 #
445 # if choice==0: # pick the single gaussian fit
446 # Amplitude0=lsq1[0][2]
447 # shift0=lsq1[0][0]
448 # width0=lsq1[0][1]
449 # p0=lsq1[0][3]
450 # Amplitude1=0.
451 # shift1=0.
452 # width1=0.
453 # p1=0.
454 # noise=lsq1[0][4]
455 # elif choice==1: # take the first one of the 2 gaussians fitted
456 # Amplitude0 = lsq2[0][2]
457 # shift0 = lsq2[0][0]
458 # width0 = lsq2[0][1]
459 # p0 = lsq2[0][3]
460 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
461 # shift1 = lsq2[0][4] # This is 0 in gg1
462 # width1 = lsq2[0][5] # This is 0 in gg1
463 # p1 = lsq2[0][7] # This is 0 in gg1
464 # noise = lsq2[0][8]
465 # else: # the second one
466 # Amplitude0 = lsq2[0][6]
467 # shift0 = lsq2[0][4]
468 # width0 = lsq2[0][5]
469 # p0 = lsq2[0][7]
470 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
471 # shift1 = lsq2[0][0] # This is 0 in gg1
472 # width1 = lsq2[0][1] # This is 0 in gg1
473 # p1 = lsq2[0][3] # This is 0 in gg1
474 # noise = lsq2[0][8]
475 #
476 # #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
477 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
478 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
479 # #print 'SPC_ch1.shape',SPC_ch1.shape
480 # #print 'SPC_ch2.shape',SPC_ch2.shape
481 # #dataOut.data_param = SPC_ch1
482 # GauSPC[0] = SPC_ch1
483 # GauSPC[1] = SPC_ch2
484
485 # #plt.gcf().clear()
486 # plt.figure(50+self.i)
487 # self.i=self.i+1
488 # #plt.subplot(121)
489 # plt.plot(self.spc,'k')#,label='spc(66)')
490 # plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
491 # #plt.plot(SPC_ch2,'r')#,label='gg2')
492 # #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
493 # #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
494 # #plt.plot(xFrec,FitGauss,'yo:',label='fit')
495 # plt.legend()
496 # plt.title('DATOS A ALTURA DE 7500 METROS')
497 # plt.show()
498 # print 'shift0', shift0
499 # print 'Amplitude0', Amplitude0
500 # print 'width0', width0
501 # print 'p0', p0
502 # print '========================'
503 # print 'shift1', shift1
504 # print 'Amplitude1', Amplitude1
505 # print 'width1', width1
506 # print 'p1', p1
507 # print 'noise', noise
508 # print 's_noise', wnoise
509
510 print '========================================================'
511 print 'total_time: ', time.time()-start_time
512
513 # re-normalizing spc and noise
514 # This part differs from gg1
515
516
517
518 ''' Parameters:
519 1. Amplitude
520 2. Shift
521 3. Width
522 4. Power
523 '''
524
525
526 ###############################################################################
527 def FitGau(self, X):
528
529 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
530 #print 'VARSSSS', ch, pnoise, noise, num_intg
531
532 #print 'HEIGHTS', self.Num_Hei
533
534 GauSPC = []
535 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
536 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
537 SPC_ch1[:] = 0#numpy.NaN
538 SPC_ch2[:] = 0#numpy.NaN
539
540
541
542 for ht in range(self.Num_Hei):
543 #print (numpy.asarray(self.spc).shape)
544
545 #print 'TTTTT', ch , ht
546 #print self.spc.shape
547
548
549 spc = numpy.asarray(self.spc)[ch,:,ht]
550
551 #############################################
552 # normalizing spc and noise
553 # This part differs from gg1
554 spc_norm_max = max(spc)
555 spc = spc / spc_norm_max
556 pnoise = pnoise / spc_norm_max
557 #############################################
558
559 fatspectra=1.0
560
561 wnoise = noise_ / spc_norm_max
562 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
563 #if wnoise>1.1*pnoise: # to be tested later
564 # wnoise=pnoise
565 noisebl=wnoise*0.9; noisebh=wnoise*1.1
566 spc=spc-wnoise
567 # print 'wnoise', noise_[0], spc_norm_max, wnoise
568 minx=numpy.argmin(spc)
569 spcs=numpy.roll(spc,-minx)
570 cum=numpy.cumsum(spcs)
571 tot_noise=wnoise * self.Num_Bin #64;
572 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
573 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
574 #snr=tot_signal/tot_noise
575 #snr=cum[-1]/tot_noise
576 snr = sum(spcs)/tot_noise
577 snrdB=10.*numpy.log10(snr)
578
579 if snrdB < SNRlimit :
580 snr = numpy.NaN
581 SPC_ch1[:,ht] = 0#numpy.NaN
582 SPC_ch1[:,ht] = 0#numpy.NaN
583 GauSPC = (SPC_ch1,SPC_ch2)
584 continue
585 #print 'snr',snrdB #, sum(spcs) , tot_noise
586
587
588
589 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
590 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
591
592 cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
593 cumlo=cummax*epsi;
594 cumhi=cummax*(1-epsi)
595 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
596
597
598 if len(powerindex) < 1:# case for powerindex 0
599 continue
600 powerlo=powerindex[0]
601 powerhi=powerindex[-1]
602 powerwidth=powerhi-powerlo
603
604 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
605 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
606 midpeak=(firstpeak+secondpeak)/2.
607 firstamp=spcs[int(firstpeak)]
608 secondamp=spcs[int(secondpeak)]
609 midamp=spcs[int(midpeak)]
610
611 x=numpy.arange( self.Num_Bin )
612 y_data=spc+wnoise
613
614 # single gaussian
615 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
616 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
617 power0=2.
618 amplitude0=midamp
619 state0=[shift0,width0,amplitude0,power0,wnoise]
620 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
621 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
622
623 chiSq1=lsq1[1];
624 jack1= self.y_jacobian1(x,lsq1[0])
625
626
627 try:
628 sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
629 except:
630 std1=32.; sigmas1=numpy.ones(5)
631 else:
632 std1=sigmas1[0]
633
634
635 if fatspectra<1.0 and powerwidth<4:
636 choice=0
637 Amplitude0=lsq1[0][2]
638 shift0=lsq1[0][0]
639 width0=lsq1[0][1]
640 p0=lsq1[0][3]
641 Amplitude1=0.
642 shift1=0.
643 width1=0.
644 p1=0.
645 noise=lsq1[0][4]
646 #return (numpy.array([shift0,width0,Amplitude0,p0]),
647 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
648
649 # two gaussians
650 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
651 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
652 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
653 width0=powerwidth/6.;
654 width1=width0
655 power0=2.;
656 power1=power0
657 amplitude0=firstamp;
658 amplitude1=secondamp
659 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
660 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
661 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
662 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
663
664 lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
665
666
667 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
668
669
670 try:
671 sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
672 except:
673 std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
674 else:
675 std2a=sigmas2[0]; std2b=sigmas2[4]
676
677
678
679 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
680
681 if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
682 if oneG:
683 choice=0
684 else:
685 w1=lsq2[0][1]; w2=lsq2[0][5]
686 a1=lsq2[0][2]; a2=lsq2[0][6]
687 p1=lsq2[0][3]; p2=lsq2[0][7]
688 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
689 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
690 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
691
692 if gp1>gp2:
693 if a1>0.7*a2:
694 choice=1
695 else:
696 choice=2
697 elif gp2>gp1:
698 if a2>0.7*a1:
699 choice=2
700 else:
701 choice=1
702 else:
703 choice=numpy.argmax([a1,a2])+1
704 #else:
705 #choice=argmin([std2a,std2b])+1
706
707 else: # with low SNR go to the most energetic peak
708 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
709
710
711 shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
712 shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
713
714 max_vel = 20
715
716 #first peak will be 0, second peak will be 1
717 if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
718 shift0=lsq2[0][0]
719 width0=lsq2[0][1]
720 Amplitude0=lsq2[0][2]
721 p0=lsq2[0][3]
722
723 shift1=lsq2[0][4]
724 width1=lsq2[0][5]
725 Amplitude1=lsq2[0][6]
726 p1=lsq2[0][7]
727 noise=lsq2[0][8]
728 else:
729 shift1=lsq2[0][0]
730 width1=lsq2[0][1]
731 Amplitude1=lsq2[0][2]
732 p1=lsq2[0][3]
733
734 shift0=lsq2[0][4]
735 width0=lsq2[0][5]
736 Amplitude0=lsq2[0][6]
737 p0=lsq2[0][7]
738 noise=lsq2[0][8]
739
740 if Amplitude0<0.1: # in case the peak is noise
741 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
742 if Amplitude1<0.1:
743 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
744
745
746 # if choice==0: # pick the single gaussian fit
747 # Amplitude0=lsq1[0][2]
748 # shift0=lsq1[0][0]
749 # width0=lsq1[0][1]
750 # p0=lsq1[0][3]
751 # Amplitude1=0.
752 # shift1=0.
753 # width1=0.
754 # p1=0.
755 # noise=lsq1[0][4]
756 # elif choice==1: # take the first one of the 2 gaussians fitted
757 # Amplitude0 = lsq2[0][2]
758 # shift0 = lsq2[0][0]
759 # width0 = lsq2[0][1]
760 # p0 = lsq2[0][3]
761 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
762 # shift1 = lsq2[0][4] # This is 0 in gg1
763 # width1 = lsq2[0][5] # This is 0 in gg1
764 # p1 = lsq2[0][7] # This is 0 in gg1
765 # noise = lsq2[0][8]
766 # else: # the second one
767 # Amplitude0 = lsq2[0][6]
768 # shift0 = lsq2[0][4]
769 # width0 = lsq2[0][5]
770 # p0 = lsq2[0][7]
771 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
772 # shift1 = lsq2[0][0] # This is 0 in gg1
773 # width1 = lsq2[0][1] # This is 0 in gg1
774 # p1 = lsq2[0][3] # This is 0 in gg1
775 # noise = lsq2[0][8]
776
777 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
778 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
779 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
780 #print 'SPC_ch1.shape',SPC_ch1.shape
781 #print 'SPC_ch2.shape',SPC_ch2.shape
782 #dataOut.data_param = SPC_ch1
783 GauSPC = (SPC_ch1,SPC_ch2)
784 #GauSPC[1] = SPC_ch2
785
786 # print 'shift0', shift0
787 # print 'Amplitude0', Amplitude0
788 # print 'width0', width0
789 # print 'p0', p0
790 # print '========================'
791 # print 'shift1', shift1
792 # print 'Amplitude1', Amplitude1
793 # print 'width1', width1
794 # print 'p1', p1
795 # print 'noise', noise
796 # print 's_noise', wnoise
797
798 return GauSPC
799
800
801 def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
802 y_model=self.y_model1(x,state)
803 s0,w0,a0,p0,n=state
804 e0=((x-s0)/w0)**2;
805
806 e0u=((x-s0-self.Num_Bin)/w0)**2;
807
808 e0d=((x-s0+self.Num_Bin)/w0)**2
809 m0=numpy.exp(-0.5*e0**(p0/2.));
810 m0u=numpy.exp(-0.5*e0u**(p0/2.));
811 m0d=numpy.exp(-0.5*e0d**(p0/2.))
812 JA=m0+m0u+m0d
813 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
814
815 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
816
817 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
818 jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
819 return jack1.T
820
821 def y_jacobian2(self,x,state):
822 y_model=self.y_model2(x,state)
823 s0,w0,a0,p0,s1,w1,a1,p1,n=state
824 e0=((x-s0)/w0)**2;
825
826 e0u=((x-s0- self.Num_Bin )/w0)**2;
827
828 e0d=((x-s0+ self.Num_Bin )/w0)**2
829 e1=((x-s1)/w1)**2;
830
831 e1u=((x-s1- self.Num_Bin )/w1)**2;
832
833 e1d=((x-s1+ self.Num_Bin )/w1)**2
834 m0=numpy.exp(-0.5*e0**(p0/2.));
835 m0u=numpy.exp(-0.5*e0u**(p0/2.));
836 m0d=numpy.exp(-0.5*e0d**(p0/2.))
837 m1=numpy.exp(-0.5*e1**(p1/2.));
838 m1u=numpy.exp(-0.5*e1u**(p1/2.));
839 m1d=numpy.exp(-0.5*e1d**(p1/2.))
840 JA=m0+m0u+m0d
841 JA1=m1+m1u+m1d
842 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
843 JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d)
844
845 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
846
847 JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)
848
849 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
850
851 JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2
852 jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
853 return jack2.T
854
855 def y_model1(self,x,state):
856 shift0,width0,amplitude0,power0,noise=state
857 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
858
859 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
860
861 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
862 return model0+model0u+model0d+noise
863
864 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
865 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
866 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
867
868 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
869
870 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
871 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
872
873 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
874
875 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
876 return model0+model0u+model0d+model1+model1u+model1d+noise
877
878 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
879
880 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
881
882 def misfit2(self,state,y_data,x,num_intg):
883 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
884
885
886 class PrecipitationProc(Operation):
887
888 '''
889 Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
890
891 Input:
892 self.dataOut.data_pre : SelfSpectra
893
894 Output:
895
896 self.dataOut.data_output : Reflectivity factor, rainfall Rate
897
898
899 Parameters affected:
900 '''
901
902
903 def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
904 tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
905
906 self.spc = dataOut.data_pre[0].copy()
907 self.Num_Hei = self.spc.shape[2]
908 self.Num_Bin = self.spc.shape[1]
909 self.Num_Chn = self.spc.shape[0]
910
911 Velrange = dataOut.abscissaList
912
913 if radar == "MIRA35C" :
914
915 Ze = self.dBZeMODE2(dataOut)
916
917 else:
918
919 self.Pt = Pt
920 self.Gt = Gt
921 self.Gr = Gr
922 self.Lambda = Lambda
923 self.aL = aL
924 self.tauW = tauW
925 self.ThetaT = ThetaT
926 self.ThetaR = ThetaR
927
928 RadarConstant = GetRadarConstant()
929 SPCmean = numpy.mean(self.spc,0)
930 ETA = numpy.zeros(self.Num_Hei)
931 Pr = numpy.sum(SPCmean,0)
932
933 #for R in range(self.Num_Hei):
934 # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
935
936 D_range = numpy.zeros(self.Num_Hei)
937 EqSec = numpy.zeros(self.Num_Hei)
938 del_V = numpy.zeros(self.Num_Hei)
939
940 for R in range(self.Num_Hei):
941 ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
942
943 h = R + Altitude #Range from ground to radar pulse altitude
944 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
945
946 D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
947 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
948
949 N_dist[R] = ETA[R] / SIGMA[R]
950
951 Ze = (ETA * Lambda**4) / (numpy.pi * Km)
952 Z = numpy.sum( N_dist * D_range**6 )
953 RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate
954
955
956 RR = (Ze/200)**(1/1.6)
957 dBRR = 10*numpy.log10(RR)
958
959 dBZe = 10*numpy.log10(Ze)
960 dataOut.data_output = Ze
961 dataOut.data_param = numpy.ones([2,self.Num_Hei])
962 dataOut.channelList = [0,1]
963 print 'channelList', dataOut.channelList
964 dataOut.data_param[0]=dBZe
965 dataOut.data_param[1]=dBRR
966 print 'RR SHAPE', dBRR.shape
967 print 'Ze SHAPE', dBZe.shape
968 print 'dataOut.data_param SHAPE', dataOut.data_param.shape
969
970
971 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
972
973 NPW = dataOut.NPW
974 COFA = dataOut.COFA
975
976 SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
977 RadarConst = dataOut.RadarConst
978 #frequency = 34.85*10**9
979
980 ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
981 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
982
983 ETA = numpy.sum(SNR,1)
984 print 'ETA' , ETA
985 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
986
987 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
988
989 for r in range(self.Num_Hei):
990
991 Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
992 #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
993
994 return Ze
995
996 def GetRadarConstant(self):
997
998 """
999 Constants:
1000
1001 Pt: Transmission Power dB
1002 Gt: Transmission Gain dB
1003 Gr: Reception Gain dB
1004 Lambda: Wavelenght m
1005 aL: Attenuation loses dB
1006 tauW: Width of transmission pulse s
1007 ThetaT: Transmission antenna bean angle rad
1008 ThetaR: Reception antenna beam angle rad
1009
1010 """
1011 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
1012 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
1013 RadarConstant = Numerator / Denominator
1014
1015 return RadarConstant
1016
1017
1018
1019 class FullSpectralAnalysis(Operation):
1020
1021 """
1022 Function that implements Full Spectral Analisys technique.
1023
1024 Input:
1025 self.dataOut.data_pre : SelfSpectra and CrossSPectra data
1026 self.dataOut.groupList : Pairlist of channels
1027 self.dataOut.ChanDist : Physical distance between receivers
1028
1029
1030 Output:
1031
1032 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
1033
1034
1035 Parameters affected: Winds, height range, SNR
1036
1037 """
1038 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
1039
1040 spc = dataOut.data_pre[0].copy()
1041 cspc = dataOut.data_pre[1].copy()
1042
1043 nChannel = spc.shape[0]
1044 nProfiles = spc.shape[1]
1045 nHeights = spc.shape[2]
1046
1047 pairsList = dataOut.groupList
1048 if dataOut.ChanDist is not None :
1049 ChanDist = dataOut.ChanDist
1050 else:
1051 ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
1052
1053 #print 'ChanDist', ChanDist
1054
1055 if dataOut.VelRange is not None:
1056 VelRange= dataOut.VelRange
1057 else:
1058 VelRange= dataOut.abscissaList
1059
1060 ySamples=numpy.ones([nChannel,nProfiles])
1061 phase=numpy.ones([nChannel,nProfiles])
1062 CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
1063 coherence=numpy.ones([nChannel,nProfiles])
1064 PhaseSlope=numpy.ones(nChannel)
1065 PhaseInter=numpy.ones(nChannel)
1066 dataSNR = dataOut.data_SNR
1067
1068
1069
1070 data = dataOut.data_pre
1071 noise = dataOut.noise
1072 print 'noise',noise
1073 #SNRdB = 10*numpy.log10(dataOut.data_SNR)
1074
1075 FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
1076 #SNRdBMean = []
1077
1078
1079 #for j in range(nHeights):
1080 # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
1081 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
1082
1083 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
1084
1085 velocityX=[]
1086 velocityY=[]
1087 velocityV=[]
1088
1089 dbSNR = 10*numpy.log10(dataSNR)
1090 dbSNR = numpy.average(dbSNR,0)
1091 for Height in range(nHeights):
1092
1093 [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
1094
1095 if abs(Vzon)<100. and abs(Vzon)> 0.:
1096 velocityX=numpy.append(velocityX, Vzon)#Vmag
1097
1098 else:
1099 print 'Vzon',Vzon
1100 velocityX=numpy.append(velocityX, numpy.NaN)
1101
1102 if abs(Vmer)<100. and abs(Vmer) > 0.:
1103 velocityY=numpy.append(velocityY, Vmer)#Vang
1104
1105 else:
1106 print 'Vmer',Vmer
1107 velocityY=numpy.append(velocityY, numpy.NaN)
1108
1109 if dbSNR[Height] > SNRlimit:
1110 velocityV=numpy.append(velocityV, FirstMoment[Height])
1111 else:
1112 velocityV=numpy.append(velocityV, numpy.NaN)
1113 #FirstMoment[Height]= numpy.NaN
1114 # if SNRdBMean[Height] <12:
1115 # FirstMoment[Height] = numpy.NaN
1116 # velocityX[Height] = numpy.NaN
1117 # velocityY[Height] = numpy.NaN
1118
1119
1120 data_output[0]=numpy.array(velocityX)
1121 data_output[1]=numpy.array(velocityY)
1122 data_output[2]=-velocityV#FirstMoment
1123
1124 print ' '
1125 #print 'FirstMoment'
1126 #print FirstMoment
1127 print 'velocityX',data_output[0]
1128 print ' '
1129 print 'velocityY',data_output[1]
1130 #print numpy.array(velocityY)
1131 print ' '
1132 #print 'SNR'
1133 #print 10*numpy.log10(dataOut.data_SNR)
1134 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1135 print ' '
1136
1137
1138 dataOut.data_output=data_output
1139 return
1140
1141
1142 def moving_average(self,x, N=2):
1143 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1144
1145 def gaus(self,xSamples,a,x0,sigma):
1146 return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
1147
1148 def Find(self,x,value):
1149 for index in range(len(x)):
1150 if x[index]==value:
1151 return index
1152
1153 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
1154
1155 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1156 phase=numpy.ones([spc.shape[0],spc.shape[1]])
1157 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
1158 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
1159 PhaseSlope=numpy.ones(spc.shape[0])
1160 PhaseInter=numpy.ones(spc.shape[0])
1161 xFrec=VelRange
1162
1163 '''Getting Eij and Nij'''
1164
1165 E01=ChanDist[0][0]
1166 N01=ChanDist[0][1]
1167
1168 E02=ChanDist[1][0]
1169 N02=ChanDist[1][1]
1170
1171 E12=ChanDist[2][0]
1172 N12=ChanDist[2][1]
1173
1174 z = spc.copy()
1175 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1176
1177 for i in range(spc.shape[0]):
1178
1179 '''****** Line of Data SPC ******'''
1180 zline=z[i,:,Height]
1181
1182 '''****** SPC is normalized ******'''
1183 FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
1184 FactNorm= FactNorm/numpy.sum(FactNorm)
1185
1186 SmoothSPC=self.moving_average(FactNorm,N=3)
1187
1188 xSamples = ar(range(len(SmoothSPC)))
1189 ySamples[i] = SmoothSPC
1190
1191 #dbSNR=10*numpy.log10(dataSNR)
1192 print ' '
1193 print ' '
1194 print ' '
1195
1196 #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
1197 print 'SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]
1198 print 'noise',noise
1199 print 'zline',zline.shape, zline[0:20]
1200 print 'FactNorm',FactNorm.shape, FactNorm[0:20]
1201 print 'FactNorm suma', numpy.sum(FactNorm)
1202
1203 for i in range(spc.shape[0]):
1204
1205 '''****** Line of Data CSPC ******'''
1206 cspcLine=cspc[i,:,Height].copy()
1207
1208 '''****** CSPC is normalized ******'''
1209 chan_index0 = pairsList[i][0]
1210 chan_index1 = pairsList[i][1]
1211 CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
1212
1213 CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
1214
1215 CSPCSamples[i] = CSPCNorm
1216 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1217
1218 coherence[i]= self.moving_average(coherence[i],N=2)
1219
1220 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1221
1222 print 'cspcLine', cspcLine.shape, cspcLine[0:20]
1223 print 'CSPCFactor', CSPCFactor#, CSPCFactor[0:20]
1224 print numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]
1225 print 'CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]
1226 print 'CSPCNorm suma', numpy.sum(CSPCNorm)
1227 print 'CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]
1228
1229 '''****** Getting fij width ******'''
1230
1231 yMean=[]
1232 yMean2=[]
1233
1234 for j in range(len(ySamples[1])):
1235 yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
1236
1237 '''******* Getting fitting Gaussian ******'''
1238 meanGauss=sum(xSamples*yMean) / len(xSamples)
1239 sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
1240
1241 print '****************************'
1242 print 'len(xSamples): ',len(xSamples)
1243 print 'yMean: ', yMean.shape, yMean[0:20]
1244 print 'ySamples', ySamples.shape, ySamples[0,0:20]
1245 print 'xSamples: ',xSamples.shape, xSamples[0:20]
1246
1247 print 'meanGauss',meanGauss
1248 print 'sigma',sigma
1249
1250 #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
1251 if dbSNR > SNRlimit :
1252 try:
1253 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
1254
1255 if numpy.amax(popt)>numpy.amax(yMean)*0.3:
1256 FitGauss=self.gaus(xSamples,*popt)
1257
1258 else:
1259 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1260 print 'Verificador: Dentro', Height
1261 except :#RuntimeError:
1262 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1263
1264
1265 else:
1266 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1267
1268 Maximun=numpy.amax(yMean)
1269 eMinus1=Maximun*numpy.exp(-1)#*0.8
1270
1271 HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
1272 HalfWidth= xFrec[HWpos]
1273 GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
1274 Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
1275
1276 #Vpos=FirstMoment[]
1277
1278 '''****** Getting Fij ******'''
1279
1280 GaussCenter=xFrec[GCpos]
1281 if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
1282 Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
1283 else:
1284 Fij=abs(GaussCenter-HalfWidth)+0.0000001
1285
1286 '''****** Getting Frecuency range of significant data ******'''
1287
1288 Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
1289
1290 if Rangpos<GCpos:
1291 Range=numpy.array([Rangpos,2*GCpos-Rangpos])
1292 elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
1293 Range=numpy.array([2*GCpos-Rangpos,Rangpos])
1294 else:
1295 Range = numpy.array([0,0])
1296
1297 print ' '
1298 print 'GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1)
1299 print 'Rangpos',Rangpos
1300 print 'RANGE: ', Range
1301 FrecRange=xFrec[Range[0]:Range[1]]
1302
1303 '''****** Getting SCPC Slope ******'''
1304
1305 for i in range(spc.shape[0]):
1306
1307 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
1308 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1309
1310 print 'FrecRange', len(FrecRange) , FrecRange
1311 print 'PhaseRange', len(PhaseRange), PhaseRange
1312 print ' '
1313 if len(FrecRange) == len(PhaseRange):
1314 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
1315 PhaseSlope[i]=slope
1316 PhaseInter[i]=intercept
1317 else:
1318 PhaseSlope[i]=0
1319 PhaseInter[i]=0
1320 else:
1321 PhaseSlope[i]=0
1322 PhaseInter[i]=0
1323
1324 '''Getting constant C'''
1325 cC=(Fij*numpy.pi)**2
1326
1327 '''****** Getting constants F and G ******'''
1328 MijEijNij=numpy.array([[E02,N02], [E12,N12]])
1329 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1330 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1331 MijResults=numpy.array([MijResult0,MijResult1])
1332 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1333
1334 '''****** Getting constants A, B and H ******'''
1335 W01=numpy.amax(coherence[0])
1336 W02=numpy.amax(coherence[1])
1337 W12=numpy.amax(coherence[2])
1338
1339 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1340 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1341 WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1342
1343 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1344
1345 WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1346 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1347
1348 VxVy=numpy.array([[cA,cH],[cH,cB]])
1349
1350 VxVyResults=numpy.array([-cF,-cG])
1351 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1352
1353 Vzon = Vy
1354 Vmer = Vx
1355 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1356 Vang=numpy.arctan2(Vmer,Vzon)
1357 Vver=xFrec[Vpos]
1358 print 'vzon y vmer', Vzon, Vmer
1359 return Vzon, Vmer, Vver, GaussCenter
1360
122 class SpectralMoments(Operation):
1361 class SpectralMoments(Operation):
123
1362
124 '''
1363 '''
@@ -167,20 +1406,21 class SpectralMoments(Operation):
167 dataOut.data_STD = data_param[:,3]
1406 dataOut.data_STD = data_param[:,3]
168 return
1407 return
169
1408
170 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
1409 def __calculateMoments(self, oldspec, oldfreq, n0,
1410 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
171
1411
172 if (nicoh is None): nicoh = 1
1412 if (nicoh == None): nicoh = 1
173 if (graph is None): graph = 0
1413 if (graph == None): graph = 0
174 if (smooth is None): smooth = 0
1414 if (smooth == None): smooth = 0
175 elif (self.smooth < 3): smooth = 0
1415 elif (self.smooth < 3): smooth = 0
176
1416
177 if (type1 is None): type1 = 0
1417 if (type1 == None): type1 = 0
178 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1418 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
179 if (snrth is None): snrth = -3
1419 if (snrth == None): snrth = -3
180 if (dc is None): dc = 0
1420 if (dc == None): dc = 0
181 if (aliasing is None): aliasing = 0
1421 if (aliasing == None): aliasing = 0
182 if (oldfd is None): oldfd = 0
1422 if (oldfd == None): oldfd = 0
183 if (wwauto is None): wwauto = 0
1423 if (wwauto == None): wwauto = 0
184
1424
185 if (n0 < 1.e-20): n0 = 1.e-20
1425 if (n0 < 1.e-20): n0 = 1.e-20
186
1426
@@ -247,8 +1487,8 class SpectralMoments(Operation):
247 num_pairs = len(pairslist)
1487 num_pairs = len(pairslist)
248
1488
249 vel = self.dataOut.abscissaList
1489 vel = self.dataOut.abscissaList
250 spectra = self.dataOut.data_pre[0]
1490 spectra = self.dataOut.data_pre
251 cspectra = self.dataOut.data_pre[1]
1491 cspectra = self.dataIn.data_cspc
252 delta_v = vel[1] - vel[0]
1492 delta_v = vel[1] - vel[0]
253
1493
254 #Calculating the power spectrum
1494 #Calculating the power spectrum
@@ -480,7 +1720,7 class SpectralFitting(Operation):
480 error1 = p0*numpy.nan
1720 error1 = p0*numpy.nan
481
1721
482 #Save
1722 #Save
483 if self.dataOut.data_param is None:
1723 if self.dataOut.data_param == None:
484 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1724 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
485 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1725 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
486
1726
@@ -526,6 +1766,9 class WindProfiler(Operation):
526
1766
527 n = None
1767 n = None
528
1768
1769 def __init__(self):
1770 Operation.__init__(self)
1771
529 def __calculateCosDir(self, elev, azim):
1772 def __calculateCosDir(self, elev, azim):
530 zen = (90 - elev)*numpy.pi/180
1773 zen = (90 - elev)*numpy.pi/180
531 azim = azim*numpy.pi/180
1774 azim = azim*numpy.pi/180
@@ -816,7 +2059,7 class WindProfiler(Operation):
816 self.__dataReady = True
2059 self.__dataReady = True
817 return
2060 return
818
2061
819 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2):
2062 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
820 '''
2063 '''
821 Function that implements winds estimation technique with detected meteors.
2064 Function that implements winds estimation technique with detected meteors.
822
2065
@@ -828,7 +2071,7 class WindProfiler(Operation):
828 '''
2071 '''
829 # print arrayMeteor.shape
2072 # print arrayMeteor.shape
830 #Settings
2073 #Settings
831 nInt = (heightMax - heightMin)/binkm
2074 nInt = (heightMax - heightMin)/2
832 # print nInt
2075 # print nInt
833 nInt = int(nInt)
2076 nInt = int(nInt)
834 # print nInt
2077 # print nInt
@@ -1130,11 +2373,6 class WindProfiler(Operation):
1130 hmax = kwargs['hmax']
2373 hmax = kwargs['hmax']
1131 else: hmax = 110
2374 else: hmax = 110
1132
2375
1133 if kwargs.has_key('BinKm'):
1134 binkm = kwargs['BinKm']
1135 else:
1136 binkm = 2
1137
1138 dataOut.outputInterval = nHours*3600
2376 dataOut.outputInterval = nHours*3600
1139
2377
1140 if self.__isConfig == False:
2378 if self.__isConfig == False:
@@ -1145,7 +2383,7 class WindProfiler(Operation):
1145
2383
1146 self.__isConfig = True
2384 self.__isConfig = True
1147
2385
1148 if self.__buffer is None:
2386 if self.__buffer == None:
1149 self.__buffer = dataOut.data_param
2387 self.__buffer = dataOut.data_param
1150 self.__firstdata = copy.copy(dataOut)
2388 self.__firstdata = copy.copy(dataOut)
1151
2389
@@ -1159,7 +2397,7 class WindProfiler(Operation):
1159
2397
1160 self.__initime += dataOut.outputInterval #to erase time offset
2398 self.__initime += dataOut.outputInterval #to erase time offset
1161
2399
1162 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm)
2400 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1163 dataOut.flagNoData = False
2401 dataOut.flagNoData = False
1164 self.__buffer = None
2402 self.__buffer = None
1165
2403
@@ -1187,7 +2425,7 class WindProfiler(Operation):
1187 else: mode = 'SA'
2425 else: mode = 'SA'
1188
2426
1189 #Borrar luego esto
2427 #Borrar luego esto
1190 if dataOut.groupList is None:
2428 if dataOut.groupList == None:
1191 dataOut.groupList = [(0,1),(0,2),(1,2)]
2429 dataOut.groupList = [(0,1),(0,2),(1,2)]
1192 groupList = dataOut.groupList
2430 groupList = dataOut.groupList
1193 C = 3e8
2431 C = 3e8
@@ -1209,7 +2447,7 class WindProfiler(Operation):
1209
2447
1210 self.__isConfig = True
2448 self.__isConfig = True
1211
2449
1212 if self.__buffer is None:
2450 if self.__buffer == None:
1213 self.__buffer = dataOut.data_param
2451 self.__buffer = dataOut.data_param
1214 self.__firstdata = copy.copy(dataOut)
2452 self.__firstdata = copy.copy(dataOut)
1215
2453
@@ -1235,6 +2473,8 class WindProfiler(Operation):
1235
2473
1236 class EWDriftsEstimation(Operation):
2474 class EWDriftsEstimation(Operation):
1237
2475
2476 def __init__(self):
2477 Operation.__init__(self)
1238
2478
1239 def __correctValues(self, heiRang, phi, velRadial, SNR):
2479 def __correctValues(self, heiRang, phi, velRadial, SNR):
1240 listPhi = phi.tolist()
2480 listPhi = phi.tolist()
@@ -1569,7 +2809,7 class SMDetection(Operation):
1569
2809
1570
2810
1571 #Getting Pairslist
2811 #Getting Pairslist
1572 if channelPositions is None:
2812 if channelPositions == None:
1573 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2813 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1574 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2814 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1575 meteorOps = SMOperations()
2815 meteorOps = SMOperations()
@@ -1700,7 +2940,7 class SMDetection(Operation):
1700 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
2940 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1701 dataOut.data_param = arrayParameters
2941 dataOut.data_param = arrayParameters
1702
2942
1703 if arrayParameters is None:
2943 if arrayParameters == None:
1704 dataOut.flagNoData = True
2944 dataOut.flagNoData = True
1705 else:
2945 else:
1706 dataOut.flagNoData = True
2946 dataOut.flagNoData = True
@@ -1894,7 +3134,7 class SMDetection(Operation):
1894
3134
1895 if (indDNthresh.size > 0):
3135 if (indDNthresh.size > 0):
1896 indEnd = indDNthresh[0] - 1
3136 indEnd = indDNthresh[0] - 1
1897 indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
3137 indInit = indUPthresh[j]
1898
3138
1899 meteor = powerAux[indInit:indEnd + 1]
3139 meteor = powerAux[indInit:indEnd + 1]
1900 indPeak = meteor.argmax() + indInit
3140 indPeak = meteor.argmax() + indInit
@@ -1949,19 +3189,19 class SMDetection(Operation):
1949 indSides = pairsarray[:,1]
3189 indSides = pairsarray[:,1]
1950
3190
1951 pairslist1 = list(pairslist)
3191 pairslist1 = list(pairslist)
1952 pairslist1.append((0,4))
3192 pairslist1.append((0,1))
1953 pairslist1.append((1,3))
3193 pairslist1.append((3,4))
1954
3194
1955 listMeteors1 = []
3195 listMeteors1 = []
1956 listPowerSeries = []
3196 listPowerSeries = []
1957 listVoltageSeries = []
3197 listVoltageSeries = []
1958 #volts has the war data
3198 #volts has the war data
1959
3199
1960 if frequency == 30.175e6:
3200 if frequency == 30e6:
1961 timeLag = 45*10**-3
3201 timeLag = 45*10**-3
1962 else:
3202 else:
1963 timeLag = 15*10**-3
3203 timeLag = 15*10**-3
1964 lag = int(numpy.ceil(timeLag/timeInterval))
3204 lag = numpy.ceil(timeLag/timeInterval)
1965
3205
1966 for i in range(len(listMeteors)):
3206 for i in range(len(listMeteors)):
1967
3207
@@ -1969,10 +3209,10 class SMDetection(Operation):
1969 meteorAux = numpy.zeros(16)
3209 meteorAux = numpy.zeros(16)
1970
3210
1971 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
3211 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1972 mHeight = int(listMeteors[i][0])
3212 mHeight = listMeteors[i][0]
1973 mStart = int(listMeteors[i][1])
3213 mStart = listMeteors[i][1]
1974 mPeak = int(listMeteors[i][2])
3214 mPeak = listMeteors[i][2]
1975 mEnd = int(listMeteors[i][3])
3215 mEnd = listMeteors[i][3]
1976
3216
1977 #get the volt data between the start and end times of the meteor
3217 #get the volt data between the start and end times of the meteor
1978 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
3218 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
@@ -2061,11 +3301,11 class SMDetection(Operation):
2061
3301
2062 threshError = 10
3302 threshError = 10
2063 #Depending if it is 30 or 50 MHz
3303 #Depending if it is 30 or 50 MHz
2064 if frequency == 30.175e6:
3304 if frequency == 30e6:
2065 timeLag = 45*10**-3
3305 timeLag = 45*10**-3
2066 else:
3306 else:
2067 timeLag = 15*10**-3
3307 timeLag = 15*10**-3
2068 lag = int(numpy.ceil(timeLag/timeInterval))
3308 lag = numpy.ceil(timeLag/timeInterval)
2069
3309
2070 listMeteors1 = []
3310 listMeteors1 = []
2071
3311
@@ -2121,14 +3361,14 class SMDetection(Operation):
2121 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
3361 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2122
3362
2123 pairslist1 = list(pairslist)
3363 pairslist1 = list(pairslist)
2124 pairslist1.append((0,4))
3364 pairslist1.append((0,1))
2125 pairslist1.append((1,3))
3365 pairslist1.append((3,4))
2126 numPairs = len(pairslist1)
3366 numPairs = len(pairslist1)
2127 #Time Lag
3367 #Time Lag
2128 timeLag = 45*10**-3
3368 timeLag = 45*10**-3
2129 c = 3e8
3369 c = 3e8
2130 lag = numpy.ceil(timeLag/timeInterval)
3370 lag = numpy.ceil(timeLag/timeInterval)
2131 freq = 30.175e6
3371 freq = 30e6
2132
3372
2133 listMeteors1 = []
3373 listMeteors1 = []
2134
3374
@@ -2149,7 +3389,7 class SMDetection(Operation):
2149 #Method 2
3389 #Method 2
2150 slopes = numpy.zeros(numPairs)
3390 slopes = numpy.zeros(numPairs)
2151 time = numpy.array([-2,-1,1,2])*timeInterval
3391 time = numpy.array([-2,-1,1,2])*timeInterval
2152 angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
3392 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
2153
3393
2154 #Correct phases
3394 #Correct phases
2155 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
3395 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
@@ -2236,7 +3476,7 class CorrectSMPhases(Operation):
2236 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
3476 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2237
3477
2238 meteorOps = SMOperations()
3478 meteorOps = SMOperations()
2239 if channelPositions is None:
3479 if channelPositions == None:
2240 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3480 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2241 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3481 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2242
3482
@@ -2413,7 +3653,7 class SMPhaseCalibration(Operation):
2413
3653
2414 self.__isConfig = True
3654 self.__isConfig = True
2415
3655
2416 if self.__buffer is None:
3656 if self.__buffer == None:
2417 self.__buffer = dataOut.data_param.copy()
3657 self.__buffer = dataOut.data_param.copy()
2418
3658
2419 else:
3659 else:
@@ -2469,7 +3709,6 class SMPhaseCalibration(Operation):
2469 phasesOff = phasesOff.reshape((1,phasesOff.size))
3709 phasesOff = phasesOff.reshape((1,phasesOff.size))
2470 dataOut.data_output = -phasesOff
3710 dataOut.data_output = -phasesOff
2471 dataOut.flagNoData = False
3711 dataOut.flagNoData = False
2472 dataOut.channelList = pairslist0
2473 self.__buffer = None
3712 self.__buffer = None
2474
3713
2475
3714
@@ -2803,3 +4042,4 class SMOperations():
2803 # error[indInvalid1] = 13
4042 # error[indInvalid1] = 13
2804 #
4043 #
2805 # return heights, error
4044 # return heights, error
4045 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,3 +1,5
1 import itertools
2
1 import numpy
3 import numpy
2
4
3 from jroproc_base import ProcessingUnit, Operation
5 from jroproc_base import ProcessingUnit, Operation
@@ -109,7 +111,14 class SpectraProc(ProcessingUnit):
109
111
110 if self.dataIn.type == "Spectra":
112 if self.dataIn.type == "Spectra":
111 self.dataOut.copy(self.dataIn)
113 self.dataOut.copy(self.dataIn)
114 <<<<<<< HEAD
112 # self.__selectPairs(pairsList)
115 # self.__selectPairs(pairsList)
116 =======
117 if not pairsList:
118 pairsList = itertools.combinations(self.dataOut.channelList, 2)
119 if self.dataOut.data_cspc is not None:
120 self.__selectPairs(pairsList)
121 >>>>>>> 8048843f4edfb980d0968f42f82054783b84e1cc
113 return True
122 return True
114
123
115 if self.dataIn.type == "Voltage":
124 if self.dataIn.type == "Voltage":
@@ -178,27 +187,21 class SpectraProc(ProcessingUnit):
178
187
179 def __selectPairs(self, pairsList):
188 def __selectPairs(self, pairsList):
180
189
181 if channelList == None:
190 if not pairsList:
182 return
191 return
183
192
184 pairsIndexListSelected = []
193 pairs = []
185
194 pairsIndex = []
186 for thisPair in pairsList:
187
195
188 if thisPair not in self.dataOut.pairsList:
196 for pair in pairsList:
197 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
189 continue
198 continue
199 pairs.append(pair)
200 pairsIndex.append(pairs.index(pair))
190
201
191 pairIndex = self.dataOut.pairsList.index(thisPair)
202 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
192
203 self.dataOut.pairsList = pairs
193 pairsIndexListSelected.append(pairIndex)
204 self.dataOut.pairsIndexList = pairsIndex
194
195 if not pairsIndexListSelected:
196 self.dataOut.data_cspc = None
197 self.dataOut.pairsList = []
198 return
199
200 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
201 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
202
205
203 return
206 return
204
207
Show file before | Show file after | Hide comments
@@ -1,6 +1,5
1 import sys
1 import sys
2 import numpy
2 import numpy
3 from profilehooks import profile
4 from scipy import interpolate
3 from scipy import interpolate
5 from schainpy import cSchain
4 from schainpy import cSchain
6 from jroproc_base import ProcessingUnit, Operation
5 from jroproc_base import ProcessingUnit, Operation
@@ -623,13 +622,6 class Decoder(Operation):
623
622
624 return self.datadecTime
623 return self.datadecTime
625
624
626 #@profile
627 def oldCorrelate(self, i, data, code_block):
628 profilesList = xrange(self.__nProfiles)
629 for j in profilesList:
630 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
631
632 #@profile
633 def __convolutionByBlockInTime(self, data):
625 def __convolutionByBlockInTime(self, data):
634
626
635 repetitions = self.__nProfiles / self.nCode
627 repetitions = self.__nProfiles / self.nCode
@@ -639,26 +631,11 class Decoder(Operation):
639 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
631 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
640 profilesList = xrange(self.__nProfiles)
632 profilesList = xrange(self.__nProfiles)
641
633
642 # def toVectorize(a,b):
643 # return numpy.correlate(a,b, mode='full')
644 # vectorized = numpy.vectorize(toVectorize, signature='(n),(m)->(k)')
645 for i in range(self.__nChannels):
634 for i in range(self.__nChannels):
646 # self.datadecTime[i,:,:] = numpy.array([numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] for j in profilesList ])
647 # def func(i, j):
648 # self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
649 # map(lambda j: func(i, j), range(self.__nProfiles))
650 #print data[i,:,:].shape
651 # self.datadecTime[i,:,:] = vectorized(data[i,:,:], code_block[:,:])[:,self.nBaud-1:]
652 for j in profilesList:
635 for j in profilesList:
653 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
636 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
654 # print data[i,:,:]
655 # print cSchain.correlateByBlock(data[i,:,:], code_block, 2)
656 # self.datadecTime[i,:,:] = cSchain.correlateByBlock(data[i,:,:], code_block, 2)
657 # print self.datadecTime[i,:,:]
658 #print self.datadecTime[i,:,:].shape
659 return self.datadecTime
637 return self.datadecTime
660
638
661
662 def __convolutionByBlockInFreq(self, data):
639 def __convolutionByBlockInFreq(self, data):
663
640
664 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
641 raise NotImplementedError, "Decoder by frequency fro Blocks not implemented"
Show file before | Show file after | Hide comments
@@ -7,7 +7,6 import json
7 import numpy
7 import numpy
8 import paho.mqtt.client as mqtt
8 import paho.mqtt.client as mqtt
9 import zmq
9 import zmq
10 from profilehooks import profile
11 import datetime
10 import datetime
12 from zmq.utils.monitor import recv_monitor_message
11 from zmq.utils.monitor import recv_monitor_message
13 from functools import wraps
12 from functools import wraps
@@ -16,6 +15,7 from multiprocessing import Process
16
15
17 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
16 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
18 from schainpy.model.data.jrodata import JROData
17 from schainpy.model.data.jrodata import JROData
18 from schainpy.utils import log
19
19
20 MAXNUMX = 100
20 MAXNUMX = 100
21 MAXNUMY = 100
21 MAXNUMY = 100
@@ -31,14 +31,13 def roundFloats(obj):
31 return round(obj, 2)
31 return round(obj, 2)
32
32
33 def decimate(z, MAXNUMY):
33 def decimate(z, MAXNUMY):
34 # dx = int(len(self.x)/self.__MAXNUMX) + 1
35
36 dy = int(len(z[0])/MAXNUMY) + 1
34 dy = int(len(z[0])/MAXNUMY) + 1
37
35
38 return z[::, ::dy]
36 return z[::, ::dy]
39
37
40 class throttle(object):
38 class throttle(object):
41 """Decorator that prevents a function from being called more than once every
39 '''
40 Decorator that prevents a function from being called more than once every
42 time period.
41 time period.
43 To create a function that cannot be called more than once a minute, but
42 To create a function that cannot be called more than once a minute, but
44 will sleep until it can be called:
43 will sleep until it can be called:
@@ -49,7 +48,7 class throttle(object):
49 for i in range(10):
48 for i in range(10):
50 foo()
49 foo()
51 print "This function has run %s times." % i
50 print "This function has run %s times." % i
52 """
51 '''
53
52
54 def __init__(self, seconds=0, minutes=0, hours=0):
53 def __init__(self, seconds=0, minutes=0, hours=0):
55 self.throttle_period = datetime.timedelta(
54 self.throttle_period = datetime.timedelta(
@@ -73,9 +72,169 class throttle(object):
73
72
74 return wrapper
73 return wrapper
75
74
75 class Data(object):
76 '''
77 Object to hold data to be plotted
78 '''
79
80 def __init__(self, plottypes, throttle_value):
81 self.plottypes = plottypes
82 self.throttle = throttle_value
83 self.ended = False
84 self.__times = []
85
86 def __str__(self):
87 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
88 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
89
90 def __len__(self):
91 return len(self.__times)
92
93 def __getitem__(self, key):
94 if key not in self.data:
95 raise KeyError(log.error('Missing key: {}'.format(key)))
96
97 if 'spc' in key:
98 ret = self.data[key]
99 else:
100 ret = numpy.array([self.data[key][x] for x in self.times])
101 if ret.ndim > 1:
102 ret = numpy.swapaxes(ret, 0, 1)
103 return ret
104
105 def setup(self):
106 '''
107 Configure object
108 '''
109
110 self.ended = False
111 self.data = {}
112 self.__times = []
113 self.__heights = []
114 self.__all_heights = set()
115 for plot in self.plottypes:
116 self.data[plot] = {}
117
118 def shape(self, key):
119 '''
120 Get the shape of the one-element data for the given key
121 '''
122
123 if len(self.data[key]):
124 if 'spc' in key:
125 return self.data[key].shape
126 return self.data[key][self.__times[0]].shape
127 return (0,)
128
129 def update(self, dataOut):
130 '''
131 Update data object with new dataOut
132 '''
133
134 tm = dataOut.utctime
135 if tm in self.__times:
136 return
137
138 self.parameters = getattr(dataOut, 'parameters', [])
139 self.pairs = dataOut.pairsList
140 self.channels = dataOut.channelList
141 self.xrange = (dataOut.getFreqRange(1)/1000. , dataOut.getAcfRange(1) , dataOut.getVelRange(1))
142 self.interval = dataOut.getTimeInterval()
143 self.__heights.append(dataOut.heightList)
144 self.__all_heights.update(dataOut.heightList)
145 self.__times.append(tm)
146
147 for plot in self.plottypes:
148 if plot == 'spc':
149 z = dataOut.data_spc/dataOut.normFactor
150 self.data[plot] = 10*numpy.log10(z)
151 if plot == 'cspc':
152 self.data[plot] = dataOut.data_cspc
153 if plot == 'noise':
154 self.data[plot][tm] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
155 if plot == 'rti':
156 self.data[plot][tm] = dataOut.getPower()
157 if plot == 'snr_db':
158 self.data['snr'][tm] = dataOut.data_SNR
159 if plot == 'snr':
160 self.data[plot][tm] = 10*numpy.log10(dataOut.data_SNR)
161 if plot == 'dop':
162 self.data[plot][tm] = 10*numpy.log10(dataOut.data_DOP)
163 if plot == 'mean':
164 self.data[plot][tm] = dataOut.data_MEAN
165 if plot == 'std':
166 self.data[plot][tm] = dataOut.data_STD
167 if plot == 'coh':
168 self.data[plot][tm] = dataOut.getCoherence()
169 if plot == 'phase':
170 self.data[plot][tm] = dataOut.getCoherence(phase=True)
171 if plot == 'output':
172 self.data[plot][tm] = dataOut.data_output
173 if plot == 'param':
174 self.data[plot][tm] = dataOut.data_param
175
176 def normalize_heights(self):
177 '''
178 Ensure same-dimension of the data for different heighList
179 '''
180
181 H = numpy.array(list(self.__all_heights))
182 H.sort()
183 for key in self.data:
184 shape = self.shape(key)[:-1] + H.shape
185 for tm, obj in self.data[key].items():
186 h = self.__heights[self.__times.index(tm)]
187 if H.size == h.size:
188 continue
189 index = numpy.where(numpy.in1d(H, h))[0]
190 dummy = numpy.zeros(shape) + numpy.nan
191 if len(shape) == 2:
192 dummy[:, index] = obj
193 else:
194 dummy[index] = obj
195 self.data[key][tm] = dummy
196
197 self.__heights = [H for tm in self.__times]
198
199 def jsonify(self, decimate=False):
200 '''
201 Convert data to json
202 '''
203
204 ret = {}
205 tm = self.times[-1]
206
207 for key, value in self.data:
208 if key in ('spc', 'cspc'):
209 ret[key] = roundFloats(self.data[key].to_list())
210 else:
211 ret[key] = roundFloats(self.data[key][tm].to_list())
212
213 ret['timestamp'] = tm
214 ret['interval'] = self.interval
215
216 @property
217 def times(self):
218 '''
219 Return the list of times of the current data
220 '''
221
222 ret = numpy.array(self.__times)
223 ret.sort()
224 return ret
225
226 @property
227 def heights(self):
228 '''
229 Return the list of heights of the current data
230 '''
231
232 return numpy.array(self.__heights[-1])
76
233
77 class PublishData(Operation):
234 class PublishData(Operation):
78 """Clase publish."""
235 '''
236 Operation to send data over zmq.
237 '''
79
238
80 def __init__(self, **kwargs):
239 def __init__(self, **kwargs):
81 """Inicio."""
240 """Inicio."""
@@ -87,11 +246,11 class PublishData(Operation):
87
246
88 def on_disconnect(self, client, userdata, rc):
247 def on_disconnect(self, client, userdata, rc):
89 if rc != 0:
248 if rc != 0:
90 print("Unexpected disconnection.")
249 log.warning('Unexpected disconnection.')
91 self.connect()
250 self.connect()
92
251
93 def connect(self):
252 def connect(self):
94 print 'trying to connect'
253 log.warning('trying to connect')
95 try:
254 try:
96 self.client.connect(
255 self.client.connect(
97 host=self.host,
256 host=self.host,
@@ -105,7 +264,7 class PublishData(Operation):
105 # retain=True
264 # retain=True
106 # )
265 # )
107 except:
266 except:
108 print "MQTT Conection error."
267 log.error('MQTT Conection error.')
109 self.client = False
268 self.client = False
110
269
111 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
270 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, verbose=True, **kwargs):
@@ -121,7 +280,6 class PublishData(Operation):
121 self.mqtt = kwargs.get('plottype', 0)
280 self.mqtt = kwargs.get('plottype', 0)
122 self.client = None
281 self.client = None
123 self.verbose = verbose
282 self.verbose = verbose
124 self.dataOut.firstdata = True
125 setup = []
283 setup = []
126 if mqtt is 1:
284 if mqtt is 1:
127 self.client = mqtt.Client(
285 self.client = mqtt.Client(
@@ -176,7 +334,6 class PublishData(Operation):
176 'type': self.plottype,
334 'type': self.plottype,
177 'yData': yData
335 'yData': yData
178 }
336 }
179 # print payload
180
337
181 elif self.plottype in ('rti', 'power'):
338 elif self.plottype in ('rti', 'power'):
182 data = getattr(self.dataOut, 'data_spc')
339 data = getattr(self.dataOut, 'data_spc')
@@ -230,15 +387,16 class PublishData(Operation):
230 'timestamp': 'None',
387 'timestamp': 'None',
231 'type': None
388 'type': None
232 }
389 }
233 # print 'Publishing data to {}'.format(self.host)
390
234 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
391 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
235
392
236 if self.zeromq is 1:
393 if self.zeromq is 1:
237 if self.verbose:
394 if self.verbose:
238 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
395 log.log(
396 '{} - {}'.format(self.dataOut.type, self.dataOut.datatime),
397 'Sending'
398 )
239 self.zmq_socket.send_pyobj(self.dataOut)
399 self.zmq_socket.send_pyobj(self.dataOut)
240 self.dataOut.firstdata = False
241
242
400
243 def run(self, dataOut, **kwargs):
401 def run(self, dataOut, **kwargs):
244 self.dataOut = dataOut
402 self.dataOut = dataOut
@@ -253,6 +411,7 class PublishData(Operation):
253 if self.zeromq is 1:
411 if self.zeromq is 1:
254 self.dataOut.finished = True
412 self.dataOut.finished = True
255 self.zmq_socket.send_pyobj(self.dataOut)
413 self.zmq_socket.send_pyobj(self.dataOut)
414 time.sleep(0.1)
256 self.zmq_socket.close()
415 self.zmq_socket.close()
257 if self.client:
416 if self.client:
258 self.client.loop_stop()
417 self.client.loop_stop()
@@ -281,7 +440,7 class ReceiverData(ProcessingUnit):
281 self.receiver = self.context.socket(zmq.PULL)
440 self.receiver = self.context.socket(zmq.PULL)
282 self.receiver.bind(self.address)
441 self.receiver.bind(self.address)
283 time.sleep(0.5)
442 time.sleep(0.5)
284 print '[Starting] ReceiverData from {}'.format(self.address)
443 log.success('ReceiverData from {}'.format(self.address))
285
444
286
445
287 def run(self):
446 def run(self):
@@ -291,8 +450,9 class ReceiverData(ProcessingUnit):
291 self.isConfig = True
450 self.isConfig = True
292
451
293 self.dataOut = self.receiver.recv_pyobj()
452 self.dataOut = self.receiver.recv_pyobj()
294 print '[Receiving] {} - {}'.format(self.dataOut.type,
453 log.log('{} - {}'.format(self.dataOut.type,
295 self.dataOut.datatime.ctime())
454 self.dataOut.datatime.ctime(),),
455 'Receiving')
296
456
297
457
298 class PlotterReceiver(ProcessingUnit, Process):
458 class PlotterReceiver(ProcessingUnit, Process):
@@ -306,7 +466,6 class PlotterReceiver(ProcessingUnit, Process):
306 self.mp = False
466 self.mp = False
307 self.isConfig = False
467 self.isConfig = False
308 self.isWebConfig = False
468 self.isWebConfig = False
309 self.plottypes = []
310 self.connections = 0
469 self.connections = 0
311 server = kwargs.get('server', 'zmq.pipe')
470 server = kwargs.get('server', 'zmq.pipe')
312 plot_server = kwargs.get('plot_server', 'zmq.web')
471 plot_server = kwargs.get('plot_server', 'zmq.web')
@@ -326,19 +485,13 class PlotterReceiver(ProcessingUnit, Process):
326 self.realtime = kwargs.get('realtime', False)
485 self.realtime = kwargs.get('realtime', False)
327 self.throttle_value = kwargs.get('throttle', 5)
486 self.throttle_value = kwargs.get('throttle', 5)
328 self.sendData = self.initThrottle(self.throttle_value)
487 self.sendData = self.initThrottle(self.throttle_value)
488 self.dates = []
329 self.setup()
489 self.setup()
330
490
331 def setup(self):
491 def setup(self):
332
492
333 self.data = {}
493 self.data = Data(self.plottypes, self.throttle_value)
334 self.data['times'] = []
335 for plottype in self.plottypes:
336 self.data[plottype] = {}
337 self.data['noise'] = {}
338 self.data['throttle'] = self.throttle_value
339 self.data['ENDED'] = False
340 self.isConfig = True
494 self.isConfig = True
341 self.data_web = {}
342
495
343 def event_monitor(self, monitor):
496 def event_monitor(self, monitor):
344
497
@@ -355,15 +508,13 class PlotterReceiver(ProcessingUnit, Process):
355 self.connections += 1
508 self.connections += 1
356 if evt['event'] == 512:
509 if evt['event'] == 512:
357 pass
510 pass
358 if self.connections == 0 and self.started is True:
359 self.ended = True
360
511
361 evt.update({'description': events[evt['event']]})
512 evt.update({'description': events[evt['event']]})
362
513
363 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
514 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
364 break
515 break
365 monitor.close()
516 monitor.close()
366 print("event monitor thread done!")
517 print('event monitor thread done!')
367
518
368 def initThrottle(self, throttle_value):
519 def initThrottle(self, throttle_value):
369
520
@@ -373,65 +524,16 class PlotterReceiver(ProcessingUnit, Process):
373
524
374 return sendDataThrottled
525 return sendDataThrottled
375
526
376
377 def send(self, data):
527 def send(self, data):
378 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
528 log.success('Sending {}'.format(data), self.name)
379 self.sender.send_pyobj(data)
529 self.sender.send_pyobj(data)
380
530
381
382 def update(self):
383 t = self.dataOut.utctime
384
385 if t in self.data['times']:
386 return
387
388 self.data['times'].append(t)
389 self.data['dataOut'] = self.dataOut
390
391 for plottype in self.plottypes:
392 if plottype == 'spc':
393 z = self.dataOut.data_spc/self.dataOut.normFactor
394 self.data[plottype] = 10*numpy.log10(z)
395 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
396 if plottype == 'cspc':
397 jcoherence = self.dataOut.data_cspc/numpy.sqrt(self.dataOut.data_spc*self.dataOut.data_spc)
398 self.data['cspc_coh'] = numpy.abs(jcoherence)
399 self.data['cspc_phase'] = numpy.arctan2(jcoherence.imag, jcoherence.real)*180/numpy.pi
400 if plottype == 'rti':
401 self.data[plottype][t] = self.dataOut.getPower()
402 if plottype == 'snr':
403 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
404 if plottype == 'dop':
405 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
406 if plottype == 'mean':
407 self.data[plottype][t] = self.dataOut.data_MEAN
408 if plottype == 'std':
409 self.data[plottype][t] = self.dataOut.data_STD
410 if plottype == 'coh':
411 self.data[plottype][t] = self.dataOut.getCoherence()
412 if plottype == 'phase':
413 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
414 if plottype == 'output':
415 self.data[plottype][t] = self.dataOut.data_output
416 if plottype == 'param':
417 self.data[plottype][t] = self.dataOut.data_param
418 if self.realtime:
419 self.data_web['timestamp'] = t
420 if plottype == 'spc':
421 self.data_web[plottype] = roundFloats(decimate(self.data[plottype]).tolist())
422 elif plottype == 'cspc':
423 self.data_web['cspc_coh'] = roundFloats(decimate(self.data['cspc_coh']).tolist())
424 self.data_web['cspc_phase'] = roundFloats(decimate(self.data['cspc_phase']).tolist())
425 elif plottype == 'noise':
426 self.data_web['noise'] = roundFloats(self.data['noise'][t].tolist())
427 else:
428 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
429 self.data_web['interval'] = self.dataOut.getTimeInterval()
430 self.data_web['type'] = plottype
431
432 def run(self):
531 def run(self):
433
532
434 print '[Starting] {} from {}'.format(self.name, self.address)
533 log.success(
534 'Starting from {}'.format(self.address),
535 self.name
536 )
435
537
436 self.context = zmq.Context()
538 self.context = zmq.Context()
437 self.receiver = self.context.socket(zmq.PULL)
539 self.receiver = self.context.socket(zmq.PULL)
@@ -448,39 +550,39 class PlotterReceiver(ProcessingUnit, Process):
448 else:
550 else:
449 self.sender.bind("ipc:///tmp/zmq.plots")
551 self.sender.bind("ipc:///tmp/zmq.plots")
450
552
451 time.sleep(3)
553 time.sleep(2)
452
554
453 t = Thread(target=self.event_monitor, args=(monitor,))
555 t = Thread(target=self.event_monitor, args=(monitor,))
454 t.start()
556 t.start()
455
557
456 while True:
558 while True:
457 self.dataOut = self.receiver.recv_pyobj()
559 dataOut = self.receiver.recv_pyobj()
458 # print '[Receiving] {} - {}'.format(self.dataOut.type,
560 dt = datetime.datetime.fromtimestamp(dataOut.utctime).date()
459 # self.dataOut.datatime.ctime())
561 sended = False
460
562 if dt not in self.dates:
461 self.update()
563 if self.data:
564 self.data.ended = True
565 self.send(self.data)
566 sended = True
567 self.data.setup()
568 self.dates.append(dt)
462
569
463 if self.dataOut.firstdata is True:
570 self.data.update(dataOut)
464 self.data['STARTED'] = True
465
571
466 if self.dataOut.finished is True:
572 if dataOut.finished is True:
467 self.send(self.data)
468 self.connections -= 1
573 self.connections -= 1
469 if self.connections == 0 and self.started:
574 if self.connections == 0 and dt in self.dates:
470 self.ended = True
575 self.data.ended = True
471 self.data['ENDED'] = True
472 self.send(self.data)
576 self.send(self.data)
473 self.setup()
577 self.data.setup()
474 self.started = False
475 else:
578 else:
476 if self.realtime:
579 if self.realtime:
477 self.send(self.data)
580 self.send(self.data)
478 self.sender_web.send_string(json.dumps(self.data_web))
581 # self.sender_web.send_string(self.data.jsonify())
479 else:
582 else:
583 if not sended:
480 self.sendData(self.send, self.data)
584 self.sendData(self.send, self.data)
481 self.started = True
482
585
483 self.data['STARTED'] = False
484 return
586 return
485
587
486 def sendToWeb(self):
588 def sendToWeb(self):
@@ -497,6 +599,6 class PlotterReceiver(ProcessingUnit, Process):
497 time.sleep(1)
599 time.sleep(1)
498 for kwargs in self.operationKwargs.values():
600 for kwargs in self.operationKwargs.values():
499 if 'plot' in kwargs:
601 if 'plot' in kwargs:
500 print '[Sending] Config data to web for {}'.format(kwargs['code'].upper())
602 log.success('[Sending] Config data to web for {}'.format(kwargs['code'].upper()))
501 sender_web_config.send_string(json.dumps(kwargs))
603 sender_web_config.send_string(json.dumps(kwargs))
502 self.isWebConfig = True No newline at end of file
604 self.isWebConfig = True
Show file before | Show file after | Hide comments
@@ -1,4 +1,4
1 """
1 '''
2 SCHAINPY - LOG
2 SCHAINPY - LOG
3 Simple helper for log standarization
3 Simple helper for log standarization
4 Usage:
4 Usage:
@@ -13,47 +13,32 SCHAINPY - LOG
13 which will look like this:
13 which will look like this:
14 [NEVER GONNA] - give you up
14 [NEVER GONNA] - give you up
15 with color red as background and white as foreground.
15 with color red as background and white as foreground.
16 """
16 '''
17 import os
17
18 import sys
19 import click
18 import click
20
19
21 def warning(message):
20 def warning(message, tag='Warning'):
22 click.echo(click.style('[WARNING] - ' + message, fg='yellow'))
21 click.echo(click.style('[{}] {}'.format(tag, message), fg='yellow'))
22 pass
23
23
24
25 def error(message, tag='Error'):
26 click.echo(click.style('[{}] {}'.format(tag, message), fg='red'))
27 pass
24
28
25 def error(message):
26 click.echo(click.style('[ERROR] - ' + message, fg='red', bg='black'))
27
29
30 def success(message, tag='Info'):
31 click.echo(click.style('[{}] {}'.format(tag, message), fg='green'))
32 pass
28
33
29 def success(message):
30 click.echo(click.style(message, fg='green'))
31
34
35 def log(message, tag='Info'):
36 click.echo('[{}] {}'.format(tag, message))
37 pass
32
38
33 def log(message, topic='LOG'):
34 click.echo('[{}] - {}'.format(topic, message))
35
39
36 def makelogger(topic, bg='reset', fg='reset'):
40 def makelogger(tag, bg='reset', fg='reset'):
37 def func(message):
41 def func(message):
38 click.echo(click.style('[{}] - '.format(topic.upper()) + message,
42 click.echo(click.style('[{}] {}'.format(tag.upper(), message),
39 bg=bg, fg=fg))
43 bg=bg, fg=fg))
40 return func
44 return func
41
42 class LoggerForFile():
43 def __init__(self, filename):
44 self.old_stdout=sys.stdout
45 cwd = os.getcwd()
46 self.log_file = open(os.path.join(cwd, filename), 'w+')
47 def write(self, text):
48 text = text.rstrip()
49 if not text:
50 return
51 self.log_file.write(text + '\n')
52 self.old_stdout.write(text + '\n')
53 def flush(self):
54 self.old_stdout.flush()
55
56 def logToFile(filename='log.log'):
57 logger = LoggerForFile(filename)
58 sys.stdout = logger
59
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now