##// END OF EJS Templates
schain-jc
RSS
data...
Julio Valdez -
r608:3cc55c179318
parent child
Show More
Show file before | Show file after | Hide comments
@@ -0,0 +1,188
1 # DIAS 19 Y 20 FEB 2014
2 # Comprobacion de Resultados DBS con SA
3
4 import os, sys
5
6 path = os.path.split(os.getcwd())[0]
7 path = os.path.split(path)[0]
8
9 sys.path.insert(0, path)
10
11 from schainpy.controller import Project
12
13 desc = "JASMET Experiment Test"
14 filename = "JASMETtest.xml"
15
16 controllerObj = Project()
17
18 controllerObj.setup(id = '191', name='test01', description=desc)
19
20 #Experimentos
21
22 #2014051 20 Feb 2014
23 # path = '/home/soporte/Data/JASMET/JASMET_30/2014106'
24 # pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1'
25
26 remotefolder = "/home/wmaster/graficos"
27 path = '/media/joscanoa/84A65E64A65E5730/soporte/Data/JASMET/JASMET_30/'
28
29 pathfile1 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/meteor')
30 pathfile2 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/wind')
31 pathfile3 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/phase')
32
33 pathfig = os.path.join(os.environ['HOME'],'Pictures/JASMET30/graph')
34
35 startTime = '00:00:00'
36 endTime = '23:59:59'
37 xmin ='17.0'
38 xmax = '34.0'
39
40 #------------------------------------------------------------------------------------------------
41 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
42 path=path,
43 startDate='2014/04/15',
44 endDate='2014/04/16',
45 startTime=startTime,
46 endTime=endTime,
47 online=0,
48 delay=5,
49 walk=1)
50
51 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
52
53
54 #--------------------------------------------------------------------------------------------------
55
56 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
57
58 opObj00 = procUnitConfObj0.addOperation(name='selectChannels')
59 opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist')
60
61 opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency')
62 # opObj01.addParameter(name='frequency', value='30.e6', format='float')
63 opObj01.addParameter(name='frequency', value='50.e6', format='float')
64
65 #opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
66 #--------------------------------------------------------------------------------------------------
67
68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj0.getId())
69 procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int')
70
71 opObj10 = procUnitConfObj1.addOperation(name='MeteorDetection')
72 # opObj10.addParameter(name='predefinedPhaseShifts', value='-89.5, 41.5, 0.0, -138.0, -85.5', format='floatlist')
73 opObj10.addParameter(name='predefinedPhaseShifts', value='0, 0, 0, 0, 0', format='floatlist')
74 opObj10.addParameter(name='cohDetection', value='0', format='bool')
75 opObj10.addParameter(name='noise_multiple', value='4', format='int')
76 opObj10.addParameter(name='SNRThresh', value='5', format='float')
77 opObj10.addParameter(name='phaseThresh', value='20', format='float')
78 opObj10.addParameter(name='azimuth', value='45', format='float')
79 opObj10.addParameter(name='hmin', value='68', format='float')
80 opObj10.addParameter(name='hmax', value='112', format='float')
81 opObj10.addParameter(name='saveAll', value='1', format='bool')
82
83 opObj12 = procUnitConfObj1.addOperation(name='HDF5Writer', optype='other')
84 opObj12.addParameter(name='path', value=pathfile1)
85 opObj12.addParameter(name='blocksPerFile', value='1000', format='int')
86 opObj12.addParameter(name='metadataList',value='type,inputUnit,heightList,paramInterval',format='list')
87 opObj12.addParameter(name='dataList',value='data_param',format='list')
88 opObj12.addParameter(name='mode',value='0',format='int')
89 #Tiene que ser de 3 dimensiones, append en lugar de aumentar una dimension
90
91 opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other')
92 opObj13.addParameter(name='id', value='1', format='int')
93 opObj13.addParameter(name='wintitle', value='Sky Map', format='str')
94 opObj13.addParameter(name='save', value='1', format='bool')
95 opObj13.addParameter(name='figpath', value=pathfig, format='str')
96 opObj13.addParameter(name='ftp', value='1', format='int')
97 opObj13.addParameter(name='exp_code', value='15', format='int')
98 opObj13.addParameter(name='sub_exp_code', value='1', format='int')
99
100
101 #--------------------------------------------------------------------------------------------------
102 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
103
104 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
105 opObj22.addParameter(name='technique', value='Meteors', format='str')
106 opObj22.addParameter(name='nHours', value='0.25', format='float')
107 opObj22.addParameter(name='SNRThresh', value='12.0', format='float')
108
109 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
110 opObj23.addParameter(name='id', value='2', format='int')
111 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
112 opObj23.addParameter(name='save', value='1', format='bool')
113 opObj23.addParameter(name='figpath', value = pathfig, format='str')
114 opObj23.addParameter(name='zmin', value='-120', format='int')
115 opObj23.addParameter(name='zmax', value='120', format='int')
116 # opObj12.addParameter(name='zmin_ver', value='-0.8', format='float')
117 # opObj12.addParameter(name='zmax_ver', value='0.8', format='float')
118 # opObj23.addParameter(name='SNRmin', value='-10', format='int')
119 # opObj23.addParameter(name='SNRmax', value='60', format='int')
120 # opObj23.addParameter(name='SNRthresh', value='0', format='float')
121 opObj23.addParameter(name='xmin', value=xmin, format='float')
122 opObj23.addParameter(name='xmax', value=xmax, format='float')
123 opObj23.addParameter(name='ftp', value='1', format='int')
124 opObj23.addParameter(name='exp_code', value='15', format='int')
125 opObj23.addParameter(name='sub_exp_code', value='1', format='int')
126
127 opObj24 = procUnitConfObj2.addOperation(name='HDF5Writer', optype='other')
128 opObj24.addParameter(name='path', value=pathfile2)
129 opObj24.addParameter(name='blocksPerFile', value='1000', format='int')
130 opObj24.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list')
131 opObj24.addParameter(name='dataList',value='data_output,utctime',format='list')
132
133 #--------------------------------------------------------------------------------------------------
134 procUnitConfObj3 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
135
136
137 opObj31 = procUnitConfObj3.addOperation(name='PhaseCalibration', optype='other')
138 opObj31.addParameter(name='nHours', value='0.25', format='float')
139 opObj31.addParameter(name='distances', value='-15, -15, 12, 12', format='intlist')
140 # opObj31.addParameter(name='distances', value='-25, -25, 20, 20', format='intlist')
141 opObj31.addParameter(name='pairs', value='(0,3),(1,2)', format='pairslist')
142 opObj31.addParameter(name='hmin', value='68', format='float')
143 opObj31.addParameter(name='hmax', value='112', format='float')
144
145 opObj32 = procUnitConfObj3.addOperation(name='PhasePlot', optype='other')
146 opObj32.addParameter(name='id', value='201', format='int')
147 opObj32.addParameter(name='wintitle', value='PhaseCalibration', format='str')
148 # opObj32.addParameter(name='timerange', value='300', format='int')
149 opObj32.addParameter(name='xmin', value=xmin, format='float')
150 opObj32.addParameter(name='xmax', value=xmax, format='float')
151 # opObj32.addParameter(name='xmin', value='19', format='float')
152 # opObj32.addParameter(name='xmax', value='20', format='float')
153 opObj32.addParameter(name='ymin', value='-180', format='float')
154 opObj32.addParameter(name='ymax', value='180', format='float')
155 opObj32.addParameter(name='figpath', value=pathfig, format='str')
156 opObj32.addParameter(name='ftp', value='1', format='int')
157 opObj32.addParameter(name='exp_code', value='15', format='int')
158 opObj32.addParameter(name='sub_exp_code', value='1', format='int')
159
160 opObj33 = procUnitConfObj3.addOperation(name='HDF5Writer', optype='other')
161 opObj33.addParameter(name='path', value=pathfile3)
162 opObj33.addParameter(name='blocksPerFile', value='1000', format='int')
163 opObj33.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list')
164 opObj33.addParameter(name='dataList',value='data_output,utctime',format='list')
165 # opObj25.addParameter(name='mode',value='1,0,0',format='intlist')
166
167 #--------------------------------------------------------------------------------------------------
168
169 procUnitConfObj4 = controllerObj.addProcUnit(name='SendToServer')
170 procUnitConfObj4.addParameter(name='server', value='jro-app.igp.gob.pe', format='str')
171 procUnitConfObj4.addParameter(name='username', value='wmaster', format='str')
172 procUnitConfObj4.addParameter(name='password', value='mst2010vhf', format='str')
173 procUnitConfObj4.addParameter(name='localfolder', value=pathfig, format='str')
174 procUnitConfObj4.addParameter(name='remotefolder', value=remotefolder, format='str')
175 procUnitConfObj4.addParameter(name='ext', value='.png', format='str')
176 procUnitConfObj4.addParameter(name='period', value=20, format='int')
177 procUnitConfObj4.addParameter(name='protocol', value='ftp', format='str')
178
179 #--------------------------------------------------------------------------------------------------
180
181 print "Escribiendo el archivo XML"
182 controllerObj.writeXml(filename)
183 print "Leyendo el archivo XML"
184 controllerObj.readXml(filename)
185
186 controllerObj.createObjects()
187 controllerObj.connectObjects()
188 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -0,0 +1,56
1 """
2 Se debe verficar que el disco de datos se encuentra montado en el sistema
3 """
4 import os, sys
5
6 path = os.path.split(os.getcwd())[0]
7 path = os.path.split(path)[0]
8
9 sys.path.insert(0, path)
10
11 from schainpy.controller import Project
12
13 desc = "Meteor Experiment Test"
14 filename = "meteor20130812.xml"
15
16 controllerObj = Project()
17 controllerObj.setup(id = '191', name='meteor_test01', description=desc)
18
19 path = '/home/dsuarez/.gvfs/data on 10.10.20.13/Jasmet50'
20
21
22 readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage',
23 path=path,
24 startDate='2014/04/16',
25 endDate='2014/04/16',
26 startTime='00:00:00',
27 endTime='23:59:59',
28 online=0,
29 walk=1)
30
31 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
32
33 procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId())
34
35
36 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
37
38
39 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
40 opObj11.addParameter(name='n', value='2', format='int')
41
42 opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other')
43 opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril')
44 opObj11.addParameter(name='blocksPerFile', value='100', format='int')
45 opObj11.addParameter(name='profilesPerBlock', value='200', format='int')
46
47
48
49 print "Escribiendo el archivo XML"
50 controllerObj.writeXml(filename)
51 print "Leyendo el archivo XML"
52 controllerObj.readXml(filename)
53
54 controllerObj.createObjects()
55 controllerObj.connectObjects()
56 controllerObj.run()
Show file before | Show file after | Hide comments
@@ -1069,13 +1069,13 class Parameters(JROData):
1069
1069
1070 data_SNR = None #Signal to Noise Ratio
1070 data_SNR = None #Signal to Noise Ratio
1071
1071
1072 heightRange = None #Heights
1072 # heightRange = None #Heights
1073
1073
1074 abscissaRange = None #Abscissa, can be velocities, lags or time
1074 abscissaList = None #Abscissa, can be velocities, lags or time
1075
1075
1076 noise = None #Noise Potency
1076 noise = None #Noise Potency
1077
1077
1078 # initUtcTime = None #Initial UTC time
1078 utctimeInit = None #Initial UTC time
1079
1079
1080 paramInterval = None #Time interval to calculate Parameters in seconds
1080 paramInterval = None #Time interval to calculate Parameters in seconds
1081
1081
@@ -1093,6 +1093,8 class Parameters(JROData):
1093
1093
1094 data_output = None #Out signal
1094 data_output = None #Out signal
1095
1095
1096
1097
1096 def __init__(self):
1098 def __init__(self):
1097 '''
1099 '''
1098 Constructor
1100 Constructor
@@ -1107,8 +1109,15 class Parameters(JROData):
1107
1109
1108 datatime = []
1110 datatime = []
1109
1111
1110 datatime.append(self.ltctime)
1112 if self.useLocalTime:
1111 datatime.append(self.ltctime + self.outputInterval - 1)
1113 time1 = self.utctimeInit - self.timeZone*60
1114 else:
1115 time1 = utctimeInit
1116
1117 # datatime.append(self.utctimeInit)
1118 # datatime.append(self.utctimeInit + self.outputInterval)
1119 datatime.append(time1)
1120 datatime.append(time1 + self.outputInterval)
1112
1121
1113 datatime = numpy.array(datatime)
1122 datatime = numpy.array(datatime)
1114
1123
Show file before | Show file after | Hide comments
@@ -84,9 +84,9 class Figure(Operation):
84 #raise ValueError, "(timerange) or (xmin & xmax) should be defined"
84 #raise ValueError, "(timerange) or (xmin & xmax) should be defined"
85
85
86 if timerange != None:
86 if timerange != None:
87 txmin = x[0] - x[0] % min(timerange/10, 10*60)
87 txmin = x[0] #- x[0] % min(timerange/10, 10*60)
88 else:
88 else:
89 txmin = x[0] - x[0] % 10*60
89 txmin = x[0] #- x[0] % 10*60
90
90
91 thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
91 thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
92 thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
92 thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
Show file before | Show file after | Hide comments
@@ -198,7 +198,7 class SkyMapPlot(Figure):
198
198
199 WIDTHPROF = None
199 WIDTHPROF = None
200 HEIGHTPROF = None
200 HEIGHTPROF = None
201 PREFIX = 'prm'
201 PREFIX = 'mmap'
202
202
203 def __init__(self):
203 def __init__(self):
204
204
@@ -213,7 +213,7 class SkyMapPlot(Figure):
213 self.HEIGHTPROF = 0
213 self.HEIGHTPROF = 0
214 self.counter_imagwr = 0
214 self.counter_imagwr = 0
215
215
216 self.PLOT_CODE = SKYMAP_CODE
216 self.PLOT_CODE = MSKYMAP_CODE
217
217
218 self.FTP_WEI = None
218 self.FTP_WEI = None
219 self.EXP_CODE = None
219 self.EXP_CODE = None
@@ -269,7 +269,7 class SkyMapPlot(Figure):
269 zmax : None
269 zmax : None
270 """
270 """
271
271
272 arrayParameters = dataOut.data_param
272 arrayParameters = dataOut.data_param[0,:]
273 error = arrayParameters[:,-1]
273 error = arrayParameters[:,-1]
274 indValid = numpy.where(error == 0)[0]
274 indValid = numpy.where(error == 0)[0]
275 finalMeteor = arrayParameters[indValid,:]
275 finalMeteor = arrayParameters[indValid,:]
@@ -325,6 +325,7 class SkyMapPlot(Figure):
325 wr_period=wr_period,
325 wr_period=wr_period,
326 thisDatetime=thisDatetime)
326 thisDatetime=thisDatetime)
327
327
328
328 class WindProfilerPlot(Figure):
329 class WindProfilerPlot(Figure):
329
330
330 __isConfig = None
331 __isConfig = None
@@ -336,7 +337,7 class WindProfilerPlot(Figure):
336
337
337 def __init__(self):
338 def __init__(self):
338
339
339 self.timerange = 2*60*60
340 self.timerange = None
340 self.__isConfig = False
341 self.__isConfig = False
341 self.__nsubplots = 1
342 self.__nsubplots = 1
342
343
@@ -391,7 +392,7 class WindProfilerPlot(Figure):
391 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
392 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
392 counter += 1
393 counter += 1
393
394
394 def run(self, dataOut, id, wintitle="", channelList=None,
395 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False',
395 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
396 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
396 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
397 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
397 timerange=None, SNRthresh = None,
398 timerange=None, SNRthresh = None,
@@ -423,11 +424,11 class WindProfilerPlot(Figure):
423 raise ValueError, "Channel %d is not in dataOut.channelList"
424 raise ValueError, "Channel %d is not in dataOut.channelList"
424 channelIndexList.append(dataOut.channelList.index(channel))
425 channelIndexList.append(dataOut.channelList.index(channel))
425
426
426 if timerange != None:
427 # if timerange != None:
427 self.timerange = timerange
428 # self.timerange = timerange
428
429 #
429 tmin = None
430 # tmin = None
430 tmax = None
431 # tmax = None
431
432
432 x = dataOut.getTimeRange1()
433 x = dataOut.getTimeRange1()
433 # y = dataOut.heightList
434 # y = dataOut.heightList
@@ -453,7 +454,7 class WindProfilerPlot(Figure):
453 z[i,ind] = numpy.nan
454 z[i,ind] = numpy.nan
454
455
455
456
456 showprofile = False
457 # showprofile = False
457 # thisDatetime = dataOut.datatime
458 # thisDatetime = dataOut.datatime
458 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
459 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
459 title = wintitle + "Wind"
460 title = wintitle + "Wind"
@@ -468,6 +469,9 class WindProfilerPlot(Figure):
468 showprofile=showprofile,
469 showprofile=showprofile,
469 show=show)
470 show=show)
470
471
472 if timerange != None:
473 self.timerange = timerange
474
471 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
475 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
472
476
473 if ymin == None: ymin = numpy.nanmin(y)
477 if ymin == None: ymin = numpy.nanmin(y)
@@ -485,6 +489,7 class WindProfilerPlot(Figure):
485 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
489 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
486 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
490 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
487
491
492
488 self.FTP_WEI = ftp_wei
493 self.FTP_WEI = ftp_wei
489 self.EXP_CODE = exp_code
494 self.EXP_CODE = exp_code
490 self.SUB_EXP_CODE = sub_exp_code
495 self.SUB_EXP_CODE = sub_exp_code
@@ -492,7 +497,7 class WindProfilerPlot(Figure):
492
497
493 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
498 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
494 self.__isConfig = True
499 self.__isConfig = True
495
500 self.figfile = figfile
496
501
497 self.setWinTitle(title)
502 self.setWinTitle(title)
498
503
@@ -620,7 +625,7 class ParametersPlot(Figure):
620 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
625 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
621 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
626 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
622 parameterIndex = None, onlyPositive = False,
627 parameterIndex = None, onlyPositive = False,
623 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None,
628 SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False,
624 DOP = True,
629 DOP = True,
625 zlabel = "", parameterName = "", parameterObject = "data_param",
630 zlabel = "", parameterName = "", parameterObject = "data_param",
626 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
631 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
@@ -695,6 +700,8 class ParametersPlot(Figure):
695 xlabel = ""
700 xlabel = ""
696 ylabel = "Range (Km)"
701 ylabel = "Range (Km)"
697
702
703 if (SNR and not onlySNR): nplots = 2*nplots
704
698 if onlyPositive:
705 if onlyPositive:
699 colormap = "jet"
706 colormap = "jet"
700 zmin = 0
707 zmin = 0
@@ -735,8 +742,22 class ParametersPlot(Figure):
735
742
736 for i in range(nchan):
743 for i in range(nchan):
737
744
745 if (SNR and not onlySNR): j = 2*i
746 else: j = i
747
738 j = nGraphsByChannel*i
748 j = nGraphsByChannel*i
739
749
750 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
751 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
752
753 if not onlySNR:
754 axes = self.axesList[j*self.__nsubplots]
755 z1 = z[i,:].reshape((1,-1))
756 axes.pcolorbuffer(x, y, z1,
757 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
758 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
759 ticksize=9, cblabel=zlabel, cbsize="1%")
760
740 if DOP:
761 if DOP:
741 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
762 title = "%s Channel %d: %s" %(parameterName, channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
742
763
@@ -751,7 +772,12 class ParametersPlot(Figure):
751
772
752 if SNR:
773 if SNR:
753 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
774 title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
775 axes = self.axesList[(j)*self.__nsubplots]
776 if not onlySNR:
777 axes = self.axesList[(j + 1)*self.__nsubplots]
778
754 axes = self.axesList[(j + nGraphsByChannel-1)]
779 axes = self.axesList[(j + nGraphsByChannel-1)]
780
755 z1 = SNRdB[i,:].reshape((1,-1))
781 z1 = SNRdB[i,:].reshape((1,-1))
756 axes.pcolorbuffer(x, y, z1,
782 axes.pcolorbuffer(x, y, z1,
757 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
783 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
@@ -1155,6 +1181,178 class EWDriftsPlot(Figure):
1155 self.__isConfig = False
1181 self.__isConfig = False
1156 self.figfile = None
1182 self.figfile = None
1157
1183
1184
1185
1186
1187 class PhasePlot(Figure):
1188
1189 __isConfig = None
1190 __nsubplots = None
1191
1192 PREFIX = 'mphase'
1193
1194 def __init__(self):
1195
1196 self.timerange = 24*60*60
1197 self.__isConfig = False
1198 self.__nsubplots = 1
1199 self.counter_imagwr = 0
1200 self.WIDTH = 600
1201 self.HEIGHT = 300
1202 self.WIDTHPROF = 120
1203 self.HEIGHTPROF = 0
1204 self.xdata = None
1205 self.ydata = None
1206
1207 self.PLOT_CODE = MPHASE_CODE
1208
1209 self.FTP_WEI = None
1210 self.EXP_CODE = None
1211 self.SUB_EXP_CODE = None
1212 self.PLOT_POS = None
1213
1214
1215 self.filename_phase = None
1216
1217 self.figfile = None
1218
1219 def getSubplots(self):
1220
1221 ncol = 1
1222 nrow = 1
1223
1224 return nrow, ncol
1225
1226 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1227
1228 self.__showprofile = showprofile
1229 self.nplots = nplots
1230
1231 ncolspan = 7
1232 colspan = 6
1233 self.__nsubplots = 2
1234
1235 self.createFigure(id = id,
1236 wintitle = wintitle,
1237 widthplot = self.WIDTH+self.WIDTHPROF,
1238 heightplot = self.HEIGHT+self.HEIGHTPROF,
1239 show=show)
1240
1241 nrow, ncol = self.getSubplots()
1242
1243 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1244
1245
1246 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1247 xmin=None, xmax=None, ymin=None, ymax=None,
1248 timerange=None,
1249 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
1250 server=None, folder=None, username=None, password=None,
1251 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1252
1253
1254 if timerange != None:
1255 self.timerange = timerange
1256
1257 tmin = None
1258 tmax = None
1259 x = dataOut.getTimeRange1()
1260 y = dataOut.getHeiRange()
1261
1262
1263 #thisDatetime = dataOut.datatime
1264 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1265 title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1266 xlabel = "Local Time"
1267 ylabel = "Phase"
1268
1269
1270 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1271 phase_beacon = dataOut.data_output
1272
1273
1274 if not self.__isConfig:
1275
1276 self.nplots = phase_beacon.size
1277
1278 self.setup(id=id,
1279 nplots=self.nplots,
1280 wintitle=wintitle,
1281 showprofile=showprofile,
1282 show=show)
1283
1284 tmin, tmax = self.getTimeLim(x, xmin, xmax)
1285 if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
1286 if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
1287
1288 self.FTP_WEI = ftp_wei
1289 self.EXP_CODE = exp_code
1290 self.SUB_EXP_CODE = sub_exp_code
1291 self.PLOT_POS = plot_pos
1292
1293 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1294 self.__isConfig = True
1295 self.figfile = figfile
1296 self.xdata = numpy.array([])
1297 self.ydata = numpy.array([])
1298
1299 #open file beacon phase
1300 path = '%s%03d' %(self.PREFIX, self.id)
1301 beacon_file = os.path.join(path,'%s.txt'%self.name)
1302 self.filename_phase = os.path.join(figpath,beacon_file)
1303 #self.save_phase(self.filename_phase)
1304
1305
1306 #store data beacon phase
1307 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1308
1309 self.setWinTitle(title)
1310
1311
1312 title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1313
1314 legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)]
1315
1316 axes = self.axesList[0]
1317
1318 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1319
1320 if len(self.ydata)==0:
1321 self.ydata = phase_beacon.reshape(-1,1)
1322 else:
1323 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1324
1325
1326 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1327 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
1328 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1329 XAxisAsTime=True, grid='both'
1330 )
1331
1332 self.draw()
1333
1334 if x[1] >= self.axesList[0].xmax:
1335 self.counter_imagwr = wr_period
1336 del self.xdata
1337 del self.ydata
1338 self.__isConfig = False
1339
1340 if self.figfile == None:
1341 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1342 self.figfile = self.getFilename(name = str_datetime)
1343
1344 if figpath != '':
1345 self.counter_imagwr += 1
1346 if (self.counter_imagwr>=wr_period):
1347 # store png plot to local folder
1348 self.saveFigure(figpath, self.figfile)
1349 # store png plot to FTP server according to RT-Web format
1350 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1351 ftp_filename = os.path.join(figpath, name)
1352 self.saveFigure(figpath, ftp_filename)
1353 self.counter_imagwr = 0
1354 self.figfile = None
1355
1158 self.save(figpath=figpath,
1356 self.save(figpath=figpath,
1159 figfile=figfile,
1357 figfile=figfile,
1160 save=save,
1358 save=save,
Show file before | Show file after | Hide comments
@@ -395,9 +395,11 def createPolar(ax, x, y,
395 # ax.set_rmax(90)
395 # ax.set_rmax(90)
396 ax.set_ylim(0,90)
396 ax.set_ylim(0,90)
397 ax.set_yticks(numpy.arange(0,90,20))
397 ax.set_yticks(numpy.arange(0,90,20))
398 ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
398 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
399 # ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11')
399 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
400 # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
400 printLabels(ax, xlabel, '', title)
401 ax.yaxis.labelpad = 230
402 printLabels(ax, xlabel, ylabel, title)
401 iplot = ax.lines[-1]
403 iplot = ax.lines[-1]
402
404
403 if '0.' in matplotlib.__version__[0:2]:
405 if '0.' in matplotlib.__version__[0:2]:
@@ -423,7 +425,7 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
423 ax = iplot.get_axes()
425 ax = iplot.get_axes()
424
426
425 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
427 # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
426 printLabels(ax, xlabel, '', title)
428 printLabels(ax, xlabel, ylabel, title)
427
429
428 set_linedata(ax, x, y, idline=0)
430 set_linedata(ax, x, y, idline=0)
429
431
Show file before | Show file after | Hide comments
@@ -18,10 +18,11 NOISE_CODE = 17
18 BEACON_CODE = 18
18 BEACON_CODE = 18
19
19
20 #USED IN jroplot_parameters.py
20 #USED IN jroplot_parameters.py
21
22 MOMENTS_CODE = 20
23 SKYMAP_CODE = 21
24 WIND_CODE = 22
21 WIND_CODE = 22
25 PARMS_CODE = 23
22 MSKYMAP_CODE = 23
26 SPECFIT_CODE = 24
23 MPHASE_CODE = 24
27 EWDRIFT_CODE = 25
24
25 MOMENTS_CODE = 25
26 PARMS_CODE = 26
27 SPECFIT_CODE = 27
28 EWDRIFT_CODE = 28
Show file before | Show file after | Hide comments
@@ -3,6 +3,7 import time
3 import os
3 import os
4 import h5py
4 import h5py
5 import re
5 import re
6 import tables
6
7
7 from schainpy.model.data.jrodata import *
8 from schainpy.model.data.jrodata import *
8 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
@@ -156,7 +157,7 class HDF5Reader(ProcessingUnit):
156 for thisPath in os.listdir(single_path):
157 for thisPath in os.listdir(single_path):
157 if not os.path.isdir(os.path.join(single_path,thisPath)):
158 if not os.path.isdir(os.path.join(single_path,thisPath)):
158 continue
159 continue
159 if not isRadarFolder(thisPath):
160 if not isDoyFolder(thisPath):
160 continue
161 continue
161
162
162 dirList.append(thisPath)
163 dirList.append(thisPath)
@@ -549,6 +550,8 class HDF5Writer(Operation):
549
550
550 metaFile = None
551 metaFile = None
551
552
553 filename = None
554
552 path = None
555 path = None
553
556
554 setFile = None
557 setFile = None
@@ -589,6 +592,8 class HDF5Writer(Operation):
589
592
590 nDims = None #Number Dimensions in each dataset
593 nDims = None #Number Dimensions in each dataset
591
594
595 nDimsForDs = None
596
592 def __init__(self):
597 def __init__(self):
593
598
594 Operation.__init__(self)
599 Operation.__init__(self)
@@ -675,16 +680,26 class HDF5Writer(Operation):
675 setFile = self.setFile
680 setFile = self.setFile
676
681
677 timeTuple = time.localtime(self.dataOut.utctime)
682 timeTuple = time.localtime(self.dataOut.utctime)
683
678 subfolder = ''
684 subfolder = ''
685 fullpath = os.path.join( path, subfolder )
679
686
687 if not( os.path.exists(fullpath) ):
688 os.mkdir(fullpath)
689 setFile = -1 #inicializo mi contador de seteo
690
691 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
680 fullpath = os.path.join( path, subfolder )
692 fullpath = os.path.join( path, subfolder )
693
681 if not( os.path.exists(fullpath) ):
694 if not( os.path.exists(fullpath) ):
682 os.mkdir(fullpath)
695 os.mkdir(fullpath)
683 setFile = -1 #inicializo mi contador de seteo
696 setFile = -1 #inicializo mi contador de seteo
697
684 else:
698 else:
685 filesList = os.listdir( fullpath )
699 filesList = os.listdir( fullpath )
686 if len( filesList ) > 0:
687 filesList = sorted( filesList, key=str.lower )
700 filesList = sorted( filesList, key=str.lower )
701 if len( filesList ) > 0:
702 filesList = [k for k in filesList if 'M' in k]
688 filen = filesList[-1]
703 filen = filesList[-1]
689 # el filename debera tener el siguiente formato
704 # el filename debera tener el siguiente formato
690 # 0 1234 567 89A BCDE (hex)
705 # 0 1234 567 89A BCDE (hex)
@@ -727,18 +742,14 class HDF5Writer(Operation):
727 setFile = self.setFile
742 setFile = self.setFile
728 mode = self.mode
743 mode = self.mode
729
744
730 if self.fp != None:
731 self.fp.close()
732
733 timeTuple = time.localtime(self.dataOut.utctime)
745 timeTuple = time.localtime(self.dataOut.utctime)
734 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
746 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
735
747
736 fullpath = os.path.join( path, subfolder )
748 fullpath = os.path.join( path, subfolder )
737 if not( os.path.exists(fullpath) ):
749
738 os.mkdir(fullpath)
750 if os.path.exists(fullpath):
739 setFile = -1 #inicializo mi contador de seteo
740 else:
741 filesList = os.listdir( fullpath )
751 filesList = os.listdir( fullpath )
752 filesList = [k for k in filesList if 'D' in k]
742 if len( filesList ) > 0:
753 if len( filesList ) > 0:
743 filesList = sorted( filesList, key=str.lower )
754 filesList = sorted( filesList, key=str.lower )
744 filen = filesList[-1]
755 filen = filesList[-1]
@@ -771,17 +782,22 class HDF5Writer(Operation):
771
782
772 ds = []
783 ds = []
773 data = []
784 data = []
785 nDimsForDs = []
774
786
775 nDatas = numpy.zeros(len(self.dataList))
787 nDatas = numpy.zeros(len(self.dataList))
776 nDims = self.arrayDim[:,0]
788 nDims = self.arrayDim[:,0]
777
789
790 nDim1 = self.arrayDim[:,2]
791 nDim0 = self.arrayDim[:,3]
792
778 for i in range(len(self.dataList)):
793 for i in range(len(self.dataList)):
779
794
780 if nDims[i]==1:
795 if nDims[i]==1:
781 ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,None) , chunks = True, dtype='S20')
796 # ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype='S20')
797 ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype=numpy.float64)
782 ds.append(ds0)
798 ds.append(ds0)
783 data.append([])
799 data.append([])
784
800 nDimsForDs.append(nDims[i])
785 else:
801 else:
786
802
787 if mode[i]==0:
803 if mode[i]==0:
@@ -800,20 +816,22 class HDF5Writer(Operation):
800 tableName = strMode + str(j)
816 tableName = strMode + str(j)
801
817
802 if nDims[i] == 3:
818 if nDims[i] == 3:
803 ds0 = grp0.create_dataset(tableName, (1,1,1) , maxshape=(None,None,None), chunks=True)
819 ds0 = grp0.create_dataset(tableName, (nDim1[i],nDim0[i],1) , data = numpy.zeros((nDim1[i],nDim0[i],1)) ,maxshape=(None,nDim0[i],None), chunks=True)
804 else:
820 else:
805 ds0 = grp0.create_dataset(tableName, (1,1) , maxshape=(None,None), chunks=True)
821 ds0 = grp0.create_dataset(tableName, (1,nDim0[i]), data = numpy.zeros((1,nDim0[i])) , maxshape=(None,nDim0[i]), chunks=True)
806
822
807 ds.append(ds0)
823 ds.append(ds0)
808 data.append([])
824 data.append([])
809
825 nDimsForDs.append(nDims[i])
810 self.nDatas = nDatas
826 self.nDatas = nDatas
811 self.nDims = nDims
827 self.nDims = nDims
812
828 self.nDimsForDs = nDimsForDs
813 #Saving variables
829 #Saving variables
814 print 'Writing the file: %s'%filename
830 print 'Writing the file: %s'%filename
831 self.filename = filename
815 self.fp = fp
832 self.fp = fp
816 self.grp = grp
833 self.grp = grp
834 self.grp.attrs.modify('nRecords', 1)
817 self.ds = ds
835 self.ds = ds
818 self.data = data
836 self.data = data
819
837
@@ -823,13 +841,65 class HDF5Writer(Operation):
823 return
841 return
824
842
825 def putData(self):
843 def putData(self):
826 self.setBlock()
844
827 self.writeBlock()
845 if not self.firsttime:
846 self.fp.flush()
847 self.fp.close()
848 self.readBlock()
828
849
829 if self.blockIndex == self.blocksPerFile:
850 if self.blockIndex == self.blocksPerFile:
851
830 self.setNextFile()
852 self.setNextFile()
853
854 self.setBlock()
855 self.writeBlock()
856
857
831 return
858 return
832
859
860 def readBlock(self):
861
862 '''
863 data Array configured
864
865
866 self.data
867 '''
868 ds = self.ds
869 #Setting HDF5 File
870 fp = h5py.File(self.filename,'r+')
871 grp = fp["Data"]
872 ind = 0
873
874 # grp.attrs['blocksPerFile'] = 0
875 for i in range(len(self.dataList)):
876
877 if self.nDims[i]==1:
878 ds0 = grp[self.dataList[i]]
879 ds[ind] = ds0
880 ind += 1
881 else:
882 if self.mode[i]==0:
883 strMode = "channel"
884 else:
885 strMode = "param"
886
887 grp0 = grp[self.dataList[i]]
888
889 for j in range(int(self.nDatas[i])):
890 tableName = strMode + str(j)
891 ds0 = grp0[tableName]
892 ds[ind] = ds0
893 ind += 1
894
895
896 self.fp = fp
897 self.grp = grp
898 self.ds = ds
899
900 return
901
902
833 def setBlock(self):
903 def setBlock(self):
834 '''
904 '''
835 data Array configured
905 data Array configured
@@ -848,11 +918,11 class HDF5Writer(Operation):
848 dataAux = getattr(self.dataOut,self.dataList[i])
918 dataAux = getattr(self.dataOut,self.dataList[i])
849
919
850 if nDims[i] == 1:
920 if nDims[i] == 1:
851 data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
921 # data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
852 if not self.firsttime:
922 data[ind] = dataAux
853 data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
923 # if not self.firsttime:
924 # data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
854 ind += 1
925 ind += 1
855
856 else:
926 else:
857 for j in range(int(nDatas[i])):
927 for j in range(int(nDatas[i])):
858 if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1
928 if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1
@@ -860,17 +930,17 class HDF5Writer(Operation):
860 else:
930 else:
861 data[ind] = dataAux[:,j,:]
931 data[ind] = dataAux[:,j,:]
862
932
863 if nDims[i] == 3:
933 # if nDims[i] == 3:
864 data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
934 # data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
865
935
866 if not self.firsttime:
936 # if not self.firsttime:
867 data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
937 # data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
868
938
869 else:
939 # else:
870 data[ind] = data[ind].reshape((1,data[ind].shape[0]))
940 # data[ind] = data[ind].reshape((1,data[ind].shape[0]))
871
941
872 if not self.firsttime:
942 # if not self.firsttime:
873 data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
943 # data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
874 ind += 1
944 ind += 1
875
945
876 self.data = data
946 self.data = data
@@ -881,13 +951,45 class HDF5Writer(Operation):
881 Saves the block in the HDF5 file
951 Saves the block in the HDF5 file
882 '''
952 '''
883 for i in range(len(self.ds)):
953 for i in range(len(self.ds)):
954 if self.firsttime:
955 # self.ds[i].resize(self.data[i].shape)
956 # self.ds[i][self.blockIndex,:] = self.data[i]
957 if type(self.data[i]) == numpy.ndarray:
958 nDims1 = len(self.ds[i].shape)
959
960 if nDims1 == 3:
961 self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1))
962
884 self.ds[i].resize(self.data[i].shape)
963 self.ds[i].resize(self.data[i].shape)
885 self.ds[i][:] = self.data[i]
964 self.ds[i][:] = self.data[i]
965 else:
966 if self.nDimsForDs[i] == 1:
967 self.ds[i].resize((self.ds[i].shape[0], self.ds[i].shape[1] + 1))
968 self.ds[i][0,-1] = self.data[i]
969 elif self.nDimsForDs[i] == 2:
970 self.ds[i].resize((self.ds[i].shape[0] + 1,self.ds[i].shape[1]))
971 self.ds[i][self.blockIndex,:] = self.data[i]
972 elif self.nDimsForDs[i] == 3:
973
974 dataShape = self.data[i].shape
975 dsShape = self.ds[i].shape
976
977 if dataShape[0]==dsShape[0]:
978 self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1))
979 self.ds[i][:,:,-1] = self.data[i]
980 else:
981 self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1],self.ds[i].shape[2]))
982 self.ds[i][dsShape[0]:,:,0] = self.data[i]
983 # self.ds[i].append(self.data[i])
984 # self.fp.flush()
985 # if not self.firsttime:
986 # self.fp.root.Data._v_attrs.nRecords = self.blockIndex
886
987
887 self.blockIndex += 1
988 # if self.firsttime:
888
989 # self.fp.close()
889 self.grp.attrs.modify('nRecords', self.blockIndex)
990 # self.readBlock2()
890
991
992 self.blockIndex += 1
891 self.firsttime = False
993 self.firsttime = False
892 return
994 return
893
995
Show file before | Show file after | Hide comments
@@ -239,13 +239,34 class ParametersProc(ProcessingUnit):
239 return moments
239 return moments
240
240
241 #------------------ Get SA Parameters --------------------------
241 #------------------ Get SA Parameters --------------------------
242
242 def GetSAParameters(self):
243 def GetSAParameters(self):
243 data = self.dataOut.data_pre
244 pairslist = self.dataOut.groupList
244 crossdata = self.dataIn.data_cspc
245 num_pairs = len(pairslist)
245 a = 1
246
247 vel = self.dataOut.abscissaList
248 spectra = self.dataOut.data_pre
249 cspectra = self.dataIn.data_cspc
250 delta_v = vel[1] - vel[0]
251
252 #Calculating the power spectrum
253 spc_pow = numpy.sum(spectra, 3)*delta_v
254 #Normalizing Spectra
255 norm_spectra = spectra/spc_pow
256 #Calculating the norm_spectra at peak
257 max_spectra = numpy.max(norm_spectra, 3)
246
258
259 #Normalizing Cross Spectra
260 norm_cspectra = numpy.zeros(cspectra.shape)
247
261
262 for i in range(num_chan):
263 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
248
264
265 max_cspectra = numpy.max(norm_cspectra,2)
266 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
267
268 for i in range(num_pairs):
269 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
249 #------------------- Get Lags ----------------------------------
270 #------------------- Get Lags ----------------------------------
250
271
251 def GetLags(self):
272 def GetLags(self):
@@ -339,8 +360,8 class ParametersProc(ProcessingUnit):
339 return phase
360 return phase
340 #------------------- Detect Meteors ------------------------------
361 #------------------- Detect Meteors ------------------------------
341
362
342 def DetectMeteors(self, hei_ref = None, tauindex = 0,
363 def MeteorDetection(self, hei_ref = None, tauindex = 0,
343 predefinedPhaseShifts = None, centerReceiverIndex = 2,
364 predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
344 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
365 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
345 noise_timeStep = 4, noise_multiple = 4,
366 noise_timeStep = 4, noise_multiple = 4,
346 multDet_timeLimit = 1, multDet_rangeLimit = 3,
367 multDet_timeLimit = 1, multDet_rangeLimit = 3,
@@ -430,14 +451,21 class ParametersProc(ProcessingUnit):
430
451
431 if predefinedPhaseShifts != None:
452 if predefinedPhaseShifts != None:
432 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
453 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
433 else:
454
455 elif beaconPhaseShifts:
434 #get hardware phase shifts using beacon signal
456 #get hardware phase shifts using beacon signal
435 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
457 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
436 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
458 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
437
459
460 else:
461 hardwarePhaseShifts = numpy.zeros(5)
462
463
438 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
464 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
439 for i in range(self.dataOut.data_pre.shape[0]):
465 for i in range(self.dataOut.data_pre.shape[0]):
440 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
466 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
467
468
441 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
469 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
442
470
443 #Remove DC
471 #Remove DC
@@ -497,29 +525,40 class ParametersProc(ProcessingUnit):
497 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
525 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
498
526
499 if len(listMeteors4) > 0:
527 if len(listMeteors4) > 0:
500 #Setting New Array
501 date = repr(self.dataOut.datatime)
502 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
503
528
504 #Calculate AOA (Error N 3, 4)
505 #JONES ET AL. 1998
506 AOAthresh = numpy.pi/8
507 error = arrayParameters[:,-1]
508 phases = -arrayMeteors4[:,9:13]
509 pairsList = []
529 pairsList = []
510 pairsList.append((0,3))
530 pairx = (0,3)
511 pairsList.append((1,2))
531 pairy = (1,2)
512 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
532 pairsList.append(pairx)
533 pairsList.append(pairy)
513
534
514 #Calculate Heights (Error N 13 and 14)
535 #Setting New Array
515 error = arrayParameters[:,-1]
536 date = repr(self.dataOut.datatime)
516 Ranges = arrayParameters[:,2]
537 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
517 zenith = arrayParameters[:,5]
538
518 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
539 meteorOps = MeteorOperations()
540 jph = numpy.array([0,0,0,0])
541 h = (hmin,hmax)
542 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
543
544 # #Calculate AOA (Error N 3, 4)
545 # #JONES ET AL. 1998
546 # error = arrayParameters[:,-1]
547 # AOAthresh = numpy.pi/8
548 # phases = -arrayParameters[:,9:13]
549 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
550 #
551 # #Calculate Heights (Error N 13 and 14)
552 # error = arrayParameters[:,-1]
553 # Ranges = arrayParameters[:,2]
554 # zenith = arrayParameters[:,5]
555 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
556 # error = arrayParameters[:,-1]
519 #********************* END OF PARAMETERS CALCULATION **************************
557 #********************* END OF PARAMETERS CALCULATION **************************
520
558
521 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
559 #***************************+ PASS DATA TO NEXT STEP **********************
522 self.dataOut.data_param = arrayParameters
560 arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
561 self.dataOut.data_param = arrayFinal
523
562
524 return
563 return
525
564
@@ -948,9 +987,8 class ParametersProc(ProcessingUnit):
948 listMeteors1 = []
987 listMeteors1 = []
949
988
950 for i in range(len(listMeteors)):
989 for i in range(len(listMeteors)):
951 meteor = listMeteors[i]
990 meteorAux = listMeteors[i]
952 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
991 if meteorAux[-1] == 0:
953 if meteor[-1] == 0:
954 mStart = listMeteors[i][1]
992 mStart = listMeteors[i][1]
955 mPeak = listMeteors[i][2]
993 mPeak = listMeteors[i][2]
956 mLag = mPeak - mStart + lag
994 mLag = mPeak - mStart + lag
@@ -1007,7 +1045,7 class ParametersProc(ProcessingUnit):
1007
1045
1008 #New arrays
1046 #New arrays
1009 arrayMeteors = numpy.array(listMeteors)
1047 arrayMeteors = numpy.array(listMeteors)
1010 arrayParameters = numpy.zeros((len(listMeteors),10))
1048 arrayParameters = numpy.zeros((len(listMeteors), 14))
1011
1049
1012 #Date inclusion
1050 #Date inclusion
1013 date = re.findall(r'\((.*?)\)', date)
1051 date = re.findall(r'\((.*?)\)', date)
@@ -1017,14 +1055,18 class ParametersProc(ProcessingUnit):
1017 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1055 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1018
1056
1019 #Meteor array
1057 #Meteor array
1020 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1058 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1021 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1059 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1022
1060
1023 #Parameters Array
1061 #Parameters Array
1024 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1062 arrayParameters[:,:2] = arrayDate #Date
1025 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1063 arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range
1064 arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error
1065 arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases
1066 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
1067
1026
1068
1027 return arrayMeteors, arrayParameters
1069 return arrayParameters
1028
1070
1029 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1071 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1030
1072
@@ -1246,7 +1288,6 class ParametersProc(ProcessingUnit):
1246 self.dataOut.data_param[i,:,h] = minp
1288 self.dataOut.data_param[i,:,h] = minp
1247 return
1289 return
1248
1290
1249
1250 def __residFunction(self, p, dp, LT, constants):
1291 def __residFunction(self, p, dp, LT, constants):
1251
1292
1252 fm = self.dataOut.library.modelFunction(p, constants)
1293 fm = self.dataOut.library.modelFunction(p, constants)
@@ -1327,7 +1368,6 class WindProfiler(Operation):
1327 return A1
1368 return A1
1328
1369
1329 def __correctValues(self, heiRang, phi, velRadial, SNR):
1370 def __correctValues(self, heiRang, phi, velRadial, SNR):
1330
1331 listPhi = phi.tolist()
1371 listPhi = phi.tolist()
1332 maxid = listPhi.index(max(listPhi))
1372 maxid = listPhi.index(max(listPhi))
1333 minid = listPhi.index(min(listPhi))
1373 minid = listPhi.index(min(listPhi))
@@ -1540,13 +1580,14 class WindProfiler(Operation):
1540
1580
1541 Parameters affected: Winds
1581 Parameters affected: Winds
1542 '''
1582 '''
1583 print arrayMeteor.shape
1543 #Settings
1584 #Settings
1544 nInt = (heightMax - heightMin)/2
1585 nInt = (heightMax - heightMin)/2
1545 winds = numpy.zeros((2,nInt))*numpy.nan
1586 winds = numpy.zeros((2,nInt))*numpy.nan
1546
1587
1547 #Filter errors
1588 #Filter errors
1548 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1589 error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
1549 finalMeteor = arrayMeteor[error,:]
1590 finalMeteor = arrayMeteor[0,error,:]
1550
1591
1551 #Meteor Histogram
1592 #Meteor Histogram
1552 finalHeights = finalMeteor[:,3]
1593 finalHeights = finalMeteor[:,3]
@@ -1592,10 +1633,10 class WindProfiler(Operation):
1592 def run(self, dataOut, technique, **kwargs):
1633 def run(self, dataOut, technique, **kwargs):
1593
1634
1594 param = dataOut.data_param
1635 param = dataOut.data_param
1595 if dataOut.abscissaList != None:
1636 # if dataOut.abscissaList != None:
1596 absc = dataOut.abscissaList[:-1]
1637 # absc = dataOut.abscissaList[:-1]
1597 noise = dataOut.noise
1638 noise = dataOut.noise
1598 heightList = dataOut.getHeiRange()
1639 heightList = dataOut.heightList
1599 SNR = dataOut.data_SNR
1640 SNR = dataOut.data_SNR
1600
1641
1601 if technique == 'DBS':
1642 if technique == 'DBS':
@@ -1616,19 +1657,14 class WindProfiler(Operation):
1616 else: correctFactor = 1
1657 else: correctFactor = 1
1617 if kwargs.has_key('channelList'):
1658 if kwargs.has_key('channelList'):
1618 channelList = kwargs['channelList']
1659 channelList = kwargs['channelList']
1660 if len(channelList) == 2:
1661 horizontalOnly = True
1619 arrayChannel = numpy.array(channelList)
1662 arrayChannel = numpy.array(channelList)
1620 param = param[arrayChannel,:,:]
1663 param = param[arrayChannel,:,:]
1621 theta_x = theta_x[arrayChannel]
1664 theta_x = theta_x[arrayChannel]
1622 theta_y = theta_y[arrayChannel]
1665 theta_y = theta_y[arrayChannel]
1623
1666
1624 velRadial0 = param[:,1,:] #Radial velocity
1667 velRadial0 = param[:,1,:] #Radial velocity
1625
1626 if velRadial0.shape[0] != theta_x.shape[0] or velRadial0.shape[0] != theta_y.shape[0]:
1627 raise ValueError, "The max number of channels is %d, and the length of cosine director is %d. Please check: dirCosX, dirCosY, elevation or azimuth arguments" %(velRadial0.shape[0], theta_x.shape[0])
1628
1629 if theta_x.shape[0] == 2:
1630 horizontalOnly = True
1631
1632 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1668 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1633 dataOut.utctimeInit = dataOut.utctime
1669 dataOut.utctimeInit = dataOut.utctime
1634 dataOut.outputInterval = dataOut.timeInterval
1670 dataOut.outputInterval = dataOut.timeInterval
@@ -1685,7 +1721,7 class WindProfiler(Operation):
1685 if self.__isConfig == False:
1721 if self.__isConfig == False:
1686 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1722 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1687 #Get Initial LTC time
1723 #Get Initial LTC time
1688 self.__initime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
1724 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1689 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1725 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1690
1726
1691 self.__isConfig = True
1727 self.__isConfig = True
@@ -1695,7 +1731,7 class WindProfiler(Operation):
1695 self.__firstdata = copy.copy(dataOut)
1731 self.__firstdata = copy.copy(dataOut)
1696
1732
1697 else:
1733 else:
1698 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1734 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1699
1735
1700 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1736 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1701
1737
@@ -1779,8 +1815,330 class EWDriftsEstimation(Operation):
1779 dataOut.outputInterval = dataOut.timeInterval
1815 dataOut.outputInterval = dataOut.timeInterval
1780 return
1816 return
1781
1817
1818 class PhaseCalibration(Operation):
1819
1820 __buffer = None
1821
1822 __initime = None
1823
1824 __dataReady = False
1825
1826 __isConfig = False
1827
1828 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
1829
1830 dataTime = currentTime + paramInterval
1831 deltaTime = dataTime - initTime
1832
1833 if deltaTime >= outputInterval or deltaTime < 0:
1834 return True
1835
1836 return False
1837
1838 def __getGammas(self, pairs, k, d, phases):
1839 gammas = numpy.zeros(2)
1840
1841 for i in range(len(pairs)):
1842
1843 pairi = pairs[i]
1844
1845 #Calculating gamma
1846 jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
1847 jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
1848
1849 #Revised distribution
1850 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
1851
1852 #Histogram
1853 nBins = 64.0
1854 rmin = -0.5*numpy.pi
1855 rmax = 0.5*numpy.pi
1856 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
1857
1858 meteorsY = phaseHisto[0]
1859 phasesX = phaseHisto[1][:-1]
1860 width = phasesX[1] - phasesX[0]
1861 phasesX += width/2
1862
1863 #Gaussian aproximation
1864 bpeak = meteorsY.argmax()
1865 peak = meteorsY.max()
1866 jmin = bpeak - 5
1867 jmax = bpeak + 5 + 1
1868
1869 if jmin<0:
1870 jmin = 0
1871 jmax = 6
1872 elif jmax > meteorsY.size:
1873 jmin = meteorsY.size - 6
1874 jmax = meteorsY.size
1875
1876 x0 = numpy.array([peak,bpeak,50])
1877 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
1878
1879 #Gammas
1880 gammas[i] = coeff[0][1]
1881 # gammas[i] = bpeak
1882
1883 return gammas
1884
1885 def __residualFunction(self, coeffs, y, t):
1886
1887 return y - self.__gauss_function(t, coeffs)
1888
1889 def __gauss_function(self, t, coeffs):
1890
1891 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
1892
1893 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
1894 meteorOps = MeteorOperations()
1895 nchan = 4
1896 pairx = pairsList[0]
1897 pairy = pairsList[1]
1898 center_xangle = 0
1899 center_yangle = 0
1900 range_angle = numpy.array([8*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
1901 ntimes = len(range_angle)
1902
1903 nstepsx = 20.0
1904 nstepsy = 20.0
1905
1906 for iz in range(ntimes):
1907 min_xangle = -range_angle[iz]/2 + center_xangle
1908 max_xangle = range_angle[iz]/2 + center_xangle
1909 min_yangle = -range_angle[iz]/2 + center_yangle
1910 max_yangle = range_angle[iz]/2 + center_yangle
1911
1912 inc_x = (max_xangle-min_xangle)/nstepsx
1913 inc_y = (max_yangle-min_yangle)/nstepsy
1914
1915 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
1916 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
1917 penalty = numpy.zeros((nstepsx,nstepsy))
1918 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
1919 jph = numpy.zeros(nchan)
1920
1921 # Iterations looking for the offset
1922 for iy in range(int(nstepsy)):
1923 for ix in range(int(nstepsx)):
1924 jph[pairy[1]] = alpha_y[iy]
1925 jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]]
1926
1927 jph[pairx[1]] = alpha_x[ix]
1928 jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]]
1929
1930 jph_array[:,ix,iy] = jph
1931
1932 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph)
1933 error = meteorsArray1[:,-1]
1934 ind1 = numpy.where(error==0)[0]
1935 penalty[ix,iy] = ind1.size
1936
1937 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
1938 phOffset = jph_array[:,i,j]
1939
1940 center_xangle = phOffset[pairx[1]]
1941 center_yangle = phOffset[pairy[1]]
1942
1943 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
1944 phOffset = phOffset*180/numpy.pi
1945 return phOffset
1946
1947
1948 def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None):
1949
1950 dataOut.flagNoData = True
1951 self.__dataReady = False
1952
1953 if nHours == None:
1954 nHours = 1
1955
1956 dataOut.outputInterval = nHours*3600
1957
1958 if self.__isConfig == False:
1959 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1960 #Get Initial LTC time
1961 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1962 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1963
1964 self.__isConfig = True
1965
1966 if self.__buffer == None:
1967 self.__buffer = dataOut.data_param.copy()
1968
1969 else:
1970 self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
1971
1972 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1973
1974 if self.__dataReady:
1975 dataOut.utctimeInit = self.__initime
1976 self.__initime += dataOut.outputInterval #to erase time offset
1977
1978 freq = dataOut.frequency
1979 c = dataOut.C #m/s
1980 lamb = c/freq
1981 k = 2*numpy.pi/lamb
1982 azimuth = 0
1983 h = (hmin, hmax)
1984 pairsList = ((0,3),(1,2))
1985
1986 meteorsArray = self.__buffer[0,:,:]
1987 error = meteorsArray[:,-1]
1988 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
1989 ind1 = numpy.where(boolError)[0]
1990 meteorsArray = meteorsArray[ind1,:]
1991 meteorsArray[:,-1] = 0
1992 phases = meteorsArray[:,9:13]
1993
1994 #Calculate Gammas
1995 gammas = self.__getGammas(pairs, k, distances, phases)
1996 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
1997 #Calculate Phases
1998 phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray)
1999 phasesOff = phasesOff.reshape((1,phasesOff.size))
2000 dataOut.data_output = -phasesOff
2001 dataOut.flagNoData = False
2002 self.__buffer = None
2003
2004
2005 return
2006
2007 class MeteorOperations():
2008
2009 def __init__(self):
2010
2011 return
2012
2013 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph):
1782
2014
2015 arrayParameters = arrayParameters0.copy()
2016 hmin = h[0]
2017 hmax = h[1]
1783
2018
2019 #Calculate AOA (Error N 3, 4)
2020 #JONES ET AL. 1998
2021 AOAthresh = numpy.pi/8
2022 error = arrayParameters[:,-1]
2023 phases = -arrayParameters[:,9:13] + jph
2024 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
2025
2026 #Calculate Heights (Error N 13 and 14)
2027 error = arrayParameters[:,-1]
2028 Ranges = arrayParameters[:,2]
2029 zenith = arrayParameters[:,5]
2030 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2031
2032 #----------------------- Get Final data ------------------------------------
2033 # error = arrayParameters[:,-1]
2034 # ind1 = numpy.where(error==0)[0]
2035 # arrayParameters = arrayParameters[ind1,:]
2036
2037 return arrayParameters
2038
2039 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2040
2041 arrayAOA = numpy.zeros((phases.shape[0],3))
2042 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2043
2044 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2045 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2046 arrayAOA[:,2] = cosDirError
2047
2048 azimuthAngle = arrayAOA[:,0]
2049 zenithAngle = arrayAOA[:,1]
2050
2051 #Setting Error
2052 #Number 3: AOA not fesible
2053 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2054 error[indInvalid] = 3
2055 #Number 4: Large difference in AOAs obtained from different antenna baselines
2056 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2057 error[indInvalid] = 4
2058 return arrayAOA, error
2059
2060 def __getDirectionCosines(self, arrayPhase, pairsList):
2061
2062 #Initializing some variables
2063 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2064 ang_aux = ang_aux.reshape(1,ang_aux.size)
2065
2066 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2067 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2068
2069
2070 for i in range(2):
2071 #First Estimation
2072 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2073 #Dealias
2074 indcsi = numpy.where(phi0_aux > numpy.pi)
2075 phi0_aux[indcsi] -= 2*numpy.pi
2076 indcsi = numpy.where(phi0_aux < -numpy.pi)
2077 phi0_aux[indcsi] += 2*numpy.pi
2078 #Direction Cosine 0
2079 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2080
2081 #Most-Accurate Second Estimation
2082 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2083 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2084 #Direction Cosine 1
2085 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2086
2087 #Searching the correct Direction Cosine
2088 cosdir0_aux = cosdir0[:,i]
2089 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2090 #Minimum Distance
2091 cosDiff = (cosdir1 - cosdir0_aux)**2
2092 indcos = cosDiff.argmin(axis = 1)
2093 #Saving Value obtained
2094 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2095
2096 return cosdir0, cosdir
2097
2098 def __calculateAOA(self, cosdir, azimuth):
2099 cosdirX = cosdir[:,0]
2100 cosdirY = cosdir[:,1]
2101
2102 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2103 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2104 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2105
2106 return angles
2107
2108 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2109
2110 Ramb = 375 #Ramb = c/(2*PRF)
2111 Re = 6371 #Earth Radius
2112 heights = numpy.zeros(Ranges.shape)
2113
2114 R_aux = numpy.array([0,1,2])*Ramb
2115 R_aux = R_aux.reshape(1,R_aux.size)
2116
2117 Ranges = Ranges.reshape(Ranges.size,1)
2118
2119 Ri = Ranges + R_aux
2120 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2121
2122 #Check if there is a height between 70 and 110 km
2123 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2124 ind_h = numpy.where(h_bool == 1)[0]
2125
2126 hCorr = hi[ind_h, :]
2127 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2128
2129 hCorr = hi[ind_hCorr]
2130 heights[ind_h] = hCorr
2131
2132 #Setting Error
2133 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2134 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2135
2136 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2137 error[indInvalid2] = 14
2138 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2139 error[indInvalid1] = 13
2140
2141 return heights, error
1784
2142
1785
2143
1786 No newline at end of file
2144
General Comments 0
You need to be logged in to leave comments. Login now