##// END OF EJS Templates
data...
Julio Valdez -
r608:3cc55c179318
parent child
Show More
@@ -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
@@ -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()
@@ -1069,13 +1069,13 class Parameters(JROData):
1069 1069
1070 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 1076 noise = None #Noise Potency
1077 1077
1078 # initUtcTime = None #Initial UTC time
1078 utctimeInit = None #Initial UTC time
1079 1079
1080 1080 paramInterval = None #Time interval to calculate Parameters in seconds
1081 1081
@@ -1093,6 +1093,8 class Parameters(JROData):
1093 1093
1094 1094 data_output = None #Out signal
1095 1095
1096
1097
1096 1098 def __init__(self):
1097 1099 '''
1098 1100 Constructor
@@ -1107,8 +1109,15 class Parameters(JROData):
1107 1109
1108 1110 datatime = []
1109 1111
1110 datatime.append(self.ltctime)
1111 datatime.append(self.ltctime + self.outputInterval - 1)
1112 if self.useLocalTime:
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 1122 datatime = numpy.array(datatime)
1114 1123
@@ -84,9 +84,9 class Figure(Operation):
84 84 #raise ValueError, "(timerange) or (xmin & xmax) should be defined"
85 85
86 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 88 else:
89 txmin = x[0] - x[0] % 10*60
89 txmin = x[0] #- x[0] % 10*60
90 90
91 91 thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
92 92 thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
@@ -198,7 +198,7 class SkyMapPlot(Figure):
198 198
199 199 WIDTHPROF = None
200 200 HEIGHTPROF = None
201 PREFIX = 'prm'
201 PREFIX = 'mmap'
202 202
203 203 def __init__(self):
204 204
@@ -213,7 +213,7 class SkyMapPlot(Figure):
213 213 self.HEIGHTPROF = 0
214 214 self.counter_imagwr = 0
215 215
216 self.PLOT_CODE = SKYMAP_CODE
216 self.PLOT_CODE = MSKYMAP_CODE
217 217
218 218 self.FTP_WEI = None
219 219 self.EXP_CODE = None
@@ -269,7 +269,7 class SkyMapPlot(Figure):
269 269 zmax : None
270 270 """
271 271
272 arrayParameters = dataOut.data_param
272 arrayParameters = dataOut.data_param[0,:]
273 273 error = arrayParameters[:,-1]
274 274 indValid = numpy.where(error == 0)[0]
275 275 finalMeteor = arrayParameters[indValid,:]
@@ -325,6 +325,7 class SkyMapPlot(Figure):
325 325 wr_period=wr_period,
326 326 thisDatetime=thisDatetime)
327 327
328
328 329 class WindProfilerPlot(Figure):
329 330
330 331 __isConfig = None
@@ -336,7 +337,7 class WindProfilerPlot(Figure):
336 337
337 338 def __init__(self):
338 339
339 self.timerange = 2*60*60
340 self.timerange = None
340 341 self.__isConfig = False
341 342 self.__nsubplots = 1
342 343
@@ -391,7 +392,7 class WindProfilerPlot(Figure):
391 392 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
392 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 396 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
396 397 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
397 398 timerange=None, SNRthresh = None,
@@ -423,11 +424,11 class WindProfilerPlot(Figure):
423 424 raise ValueError, "Channel %d is not in dataOut.channelList"
424 425 channelIndexList.append(dataOut.channelList.index(channel))
425 426
426 if timerange != None:
427 self.timerange = timerange
428
429 tmin = None
430 tmax = None
427 # if timerange != None:
428 # self.timerange = timerange
429 #
430 # tmin = None
431 # tmax = None
431 432
432 433 x = dataOut.getTimeRange1()
433 434 # y = dataOut.heightList
@@ -453,7 +454,7 class WindProfilerPlot(Figure):
453 454 z[i,ind] = numpy.nan
454 455
455 456
456 showprofile = False
457 # showprofile = False
457 458 # thisDatetime = dataOut.datatime
458 459 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
459 460 title = wintitle + "Wind"
@@ -468,6 +469,9 class WindProfilerPlot(Figure):
468 469 showprofile=showprofile,
469 470 show=show)
470 471
472 if timerange != None:
473 self.timerange = timerange
474
471 475 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
472 476
473 477 if ymin == None: ymin = numpy.nanmin(y)
@@ -485,6 +489,7 class WindProfilerPlot(Figure):
485 489 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
486 490 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
487 491
492
488 493 self.FTP_WEI = ftp_wei
489 494 self.EXP_CODE = exp_code
490 495 self.SUB_EXP_CODE = sub_exp_code
@@ -492,7 +497,7 class WindProfilerPlot(Figure):
492 497
493 498 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
494 499 self.__isConfig = True
495
500 self.figfile = figfile
496 501
497 502 self.setWinTitle(title)
498 503
@@ -620,7 +625,7 class ParametersPlot(Figure):
620 625 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
621 626 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None,
622 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 629 DOP = True,
625 630 zlabel = "", parameterName = "", parameterObject = "data_param",
626 631 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
@@ -695,6 +700,8 class ParametersPlot(Figure):
695 700 xlabel = ""
696 701 ylabel = "Range (Km)"
697 702
703 if (SNR and not onlySNR): nplots = 2*nplots
704
698 705 if onlyPositive:
699 706 colormap = "jet"
700 707 zmin = 0
@@ -735,8 +742,22 class ParametersPlot(Figure):
735 742
736 743 for i in range(nchan):
737 744
745 if (SNR and not onlySNR): j = 2*i
746 else: j = i
747
738 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 761 if DOP:
741 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 773 if SNR:
753 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 779 axes = self.axesList[(j + nGraphsByChannel-1)]
780
755 781 z1 = SNRdB[i,:].reshape((1,-1))
756 782 axes.pcolorbuffer(x, y, z1,
757 783 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
@@ -1155,6 +1181,178 class EWDriftsPlot(Figure):
1155 1181 self.__isConfig = False
1156 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 1356 self.save(figpath=figpath,
1159 1357 figfile=figfile,
1160 1358 save=save,
@@ -395,9 +395,11 def createPolar(ax, x, y,
395 395 # ax.set_rmax(90)
396 396 ax.set_ylim(0,90)
397 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 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 403 iplot = ax.lines[-1]
402 404
403 405 if '0.' in matplotlib.__version__[0:2]:
@@ -423,7 +425,7 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
423 425 ax = iplot.get_axes()
424 426
425 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 430 set_linedata(ax, x, y, idline=0)
429 431
@@ -18,10 +18,11 NOISE_CODE = 17
18 18 BEACON_CODE = 18
19 19
20 20 #USED IN jroplot_parameters.py
21
22 MOMENTS_CODE = 20
23 SKYMAP_CODE = 21
24 21 WIND_CODE = 22
25 PARMS_CODE = 23
26 SPECFIT_CODE = 24
27 EWDRIFT_CODE = 25
22 MSKYMAP_CODE = 23
23 MPHASE_CODE = 24
24
25 MOMENTS_CODE = 25
26 PARMS_CODE = 26
27 SPECFIT_CODE = 27
28 EWDRIFT_CODE = 28
@@ -3,6 +3,7 import time
3 3 import os
4 4 import h5py
5 5 import re
6 import tables
6 7
7 8 from schainpy.model.data.jrodata import *
8 9 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
@@ -156,7 +157,7 class HDF5Reader(ProcessingUnit):
156 157 for thisPath in os.listdir(single_path):
157 158 if not os.path.isdir(os.path.join(single_path,thisPath)):
158 159 continue
159 if not isRadarFolder(thisPath):
160 if not isDoyFolder(thisPath):
160 161 continue
161 162
162 163 dirList.append(thisPath)
@@ -549,6 +550,8 class HDF5Writer(Operation):
549 550
550 551 metaFile = None
551 552
553 filename = None
554
552 555 path = None
553 556
554 557 setFile = None
@@ -589,6 +592,8 class HDF5Writer(Operation):
589 592
590 593 nDims = None #Number Dimensions in each dataset
591 594
595 nDimsForDs = None
596
592 597 def __init__(self):
593 598
594 599 Operation.__init__(self)
@@ -675,16 +680,26 class HDF5Writer(Operation):
675 680 setFile = self.setFile
676 681
677 682 timeTuple = time.localtime(self.dataOut.utctime)
683
678 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 692 fullpath = os.path.join( path, subfolder )
693
681 694 if not( os.path.exists(fullpath) ):
682 695 os.mkdir(fullpath)
683 696 setFile = -1 #inicializo mi contador de seteo
697
684 698 else:
685 699 filesList = os.listdir( fullpath )
686 if len( filesList ) > 0:
687 700 filesList = sorted( filesList, key=str.lower )
701 if len( filesList ) > 0:
702 filesList = [k for k in filesList if 'M' in k]
688 703 filen = filesList[-1]
689 704 # el filename debera tener el siguiente formato
690 705 # 0 1234 567 89A BCDE (hex)
@@ -727,18 +742,14 class HDF5Writer(Operation):
727 742 setFile = self.setFile
728 743 mode = self.mode
729 744
730 if self.fp != None:
731 self.fp.close()
732
733 745 timeTuple = time.localtime(self.dataOut.utctime)
734 746 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
735 747
736 748 fullpath = os.path.join( path, subfolder )
737 if not( os.path.exists(fullpath) ):
738 os.mkdir(fullpath)
739 setFile = -1 #inicializo mi contador de seteo
740 else:
749
750 if os.path.exists(fullpath):
741 751 filesList = os.listdir( fullpath )
752 filesList = [k for k in filesList if 'D' in k]
742 753 if len( filesList ) > 0:
743 754 filesList = sorted( filesList, key=str.lower )
744 755 filen = filesList[-1]
@@ -771,17 +782,22 class HDF5Writer(Operation):
771 782
772 783 ds = []
773 784 data = []
785 nDimsForDs = []
774 786
775 787 nDatas = numpy.zeros(len(self.dataList))
776 788 nDims = self.arrayDim[:,0]
777 789
790 nDim1 = self.arrayDim[:,2]
791 nDim0 = self.arrayDim[:,3]
792
778 793 for i in range(len(self.dataList)):
779 794
780 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 798 ds.append(ds0)
783 799 data.append([])
784
800 nDimsForDs.append(nDims[i])
785 801 else:
786 802
787 803 if mode[i]==0:
@@ -800,20 +816,22 class HDF5Writer(Operation):
800 816 tableName = strMode + str(j)
801 817
802 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 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 823 ds.append(ds0)
808 824 data.append([])
809
825 nDimsForDs.append(nDims[i])
810 826 self.nDatas = nDatas
811 827 self.nDims = nDims
812
828 self.nDimsForDs = nDimsForDs
813 829 #Saving variables
814 830 print 'Writing the file: %s'%filename
831 self.filename = filename
815 832 self.fp = fp
816 833 self.grp = grp
834 self.grp.attrs.modify('nRecords', 1)
817 835 self.ds = ds
818 836 self.data = data
819 837
@@ -823,13 +841,65 class HDF5Writer(Operation):
823 841 return
824 842
825 843 def putData(self):
826 self.setBlock()
827 self.writeBlock()
844
845 if not self.firsttime:
846 self.fp.flush()
847 self.fp.close()
848 self.readBlock()
828 849
829 850 if self.blockIndex == self.blocksPerFile:
851
830 852 self.setNextFile()
853
854 self.setBlock()
855 self.writeBlock()
856
857
831 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 903 def setBlock(self):
834 904 '''
835 905 data Array configured
@@ -848,11 +918,11 class HDF5Writer(Operation):
848 918 dataAux = getattr(self.dataOut,self.dataList[i])
849 919
850 920 if nDims[i] == 1:
851 data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
852 if not self.firsttime:
853 data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
921 # data[ind] = numpy.array([str(dataAux)]).reshape((1,1))
922 data[ind] = dataAux
923 # if not self.firsttime:
924 # data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind]))
854 925 ind += 1
855
856 926 else:
857 927 for j in range(int(nDatas[i])):
858 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 930 else:
861 931 data[ind] = dataAux[:,j,:]
862 932
863 if nDims[i] == 3:
864 data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
933 # if nDims[i] == 3:
934 # data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1))
865 935
866 if not self.firsttime:
867 data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
936 # if not self.firsttime:
937 # data[ind] = numpy.dstack((self.ds[ind][:], data[ind]))
868 938
869 else:
870 data[ind] = data[ind].reshape((1,data[ind].shape[0]))
939 # else:
940 # data[ind] = data[ind].reshape((1,data[ind].shape[0]))
871 941
872 if not self.firsttime:
873 data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
942 # if not self.firsttime:
943 # data[ind] = numpy.vstack((self.ds[ind][:], data[ind]))
874 944 ind += 1
875 945
876 946 self.data = data
@@ -881,13 +951,45 class HDF5Writer(Operation):
881 951 Saves the block in the HDF5 file
882 952 '''
883 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 963 self.ds[i].resize(self.data[i].shape)
885 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
888
889 self.grp.attrs.modify('nRecords', self.blockIndex)
988 # if self.firsttime:
989 # self.fp.close()
990 # self.readBlock2()
890 991
992 self.blockIndex += 1
891 993 self.firsttime = False
892 994 return
893 995
@@ -239,13 +239,34 class ParametersProc(ProcessingUnit):
239 239 return moments
240 240
241 241 #------------------ Get SA Parameters --------------------------
242
242 243 def GetSAParameters(self):
243 data = self.dataOut.data_pre
244 crossdata = self.dataIn.data_cspc
245 a = 1
244 pairslist = self.dataOut.groupList
245 num_pairs = len(pairslist)
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 270 #------------------- Get Lags ----------------------------------
250 271
251 272 def GetLags(self):
@@ -339,8 +360,8 class ParametersProc(ProcessingUnit):
339 360 return phase
340 361 #------------------- Detect Meteors ------------------------------
341 362
342 def DetectMeteors(self, hei_ref = None, tauindex = 0,
343 predefinedPhaseShifts = None, centerReceiverIndex = 2,
363 def MeteorDetection(self, hei_ref = None, tauindex = 0,
364 predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
344 365 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
345 366 noise_timeStep = 4, noise_multiple = 4,
346 367 multDet_timeLimit = 1, multDet_rangeLimit = 3,
@@ -430,14 +451,21 class ParametersProc(ProcessingUnit):
430 451
431 452 if predefinedPhaseShifts != None:
432 453 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
433 else:
454
455 elif beaconPhaseShifts:
434 456 #get hardware phase shifts using beacon signal
435 457 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
436 458 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
437 459
460 else:
461 hardwarePhaseShifts = numpy.zeros(5)
462
463
438 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 465 for i in range(self.dataOut.data_pre.shape[0]):
440 466 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
467
468
441 469 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
442 470
443 471 #Remove DC
@@ -497,29 +525,40 class ParametersProc(ProcessingUnit):
497 525 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
498 526
499 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 529 pairsList = []
510 pairsList.append((0,3))
511 pairsList.append((1,2))
512 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
530 pairx = (0,3)
531 pairy = (1,2)
532 pairsList.append(pairx)
533 pairsList.append(pairy)
513 534
514 #Calculate Heights (Error N 13 and 14)
515 error = arrayParameters[:,-1]
516 Ranges = arrayParameters[:,2]
517 zenith = arrayParameters[:,5]
518 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
535 #Setting New Array
536 date = repr(self.dataOut.datatime)
537 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
538
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 557 #********************* END OF PARAMETERS CALCULATION **************************
520 558
521 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
522 self.dataOut.data_param = arrayParameters
559 #***************************+ PASS DATA TO NEXT STEP **********************
560 arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
561 self.dataOut.data_param = arrayFinal
523 562
524 563 return
525 564
@@ -948,9 +987,8 class ParametersProc(ProcessingUnit):
948 987 listMeteors1 = []
949 988
950 989 for i in range(len(listMeteors)):
951 meteor = listMeteors[i]
952 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
953 if meteor[-1] == 0:
990 meteorAux = listMeteors[i]
991 if meteorAux[-1] == 0:
954 992 mStart = listMeteors[i][1]
955 993 mPeak = listMeteors[i][2]
956 994 mLag = mPeak - mStart + lag
@@ -1007,7 +1045,7 class ParametersProc(ProcessingUnit):
1007 1045
1008 1046 #New arrays
1009 1047 arrayMeteors = numpy.array(listMeteors)
1010 arrayParameters = numpy.zeros((len(listMeteors),10))
1048 arrayParameters = numpy.zeros((len(listMeteors), 14))
1011 1049
1012 1050 #Date inclusion
1013 1051 date = re.findall(r'\((.*?)\)', date)
@@ -1017,14 +1055,18 class ParametersProc(ProcessingUnit):
1017 1055 arrayDate = numpy.tile(date, (len(listMeteors), 1))
1018 1056
1019 1057 #Meteor array
1020 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1021 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1058 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
1059 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
1022 1060
1023 1061 #Parameters Array
1024 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1025 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1062 arrayParameters[:,:2] = arrayDate #Date
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 1071 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1030 1072
@@ -1246,7 +1288,6 class ParametersProc(ProcessingUnit):
1246 1288 self.dataOut.data_param[i,:,h] = minp
1247 1289 return
1248 1290
1249
1250 1291 def __residFunction(self, p, dp, LT, constants):
1251 1292
1252 1293 fm = self.dataOut.library.modelFunction(p, constants)
@@ -1327,7 +1368,6 class WindProfiler(Operation):
1327 1368 return A1
1328 1369
1329 1370 def __correctValues(self, heiRang, phi, velRadial, SNR):
1330
1331 1371 listPhi = phi.tolist()
1332 1372 maxid = listPhi.index(max(listPhi))
1333 1373 minid = listPhi.index(min(listPhi))
@@ -1540,13 +1580,14 class WindProfiler(Operation):
1540 1580
1541 1581 Parameters affected: Winds
1542 1582 '''
1583 print arrayMeteor.shape
1543 1584 #Settings
1544 1585 nInt = (heightMax - heightMin)/2
1545 1586 winds = numpy.zeros((2,nInt))*numpy.nan
1546 1587
1547 1588 #Filter errors
1548 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1549 finalMeteor = arrayMeteor[error,:]
1589 error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
1590 finalMeteor = arrayMeteor[0,error,:]
1550 1591
1551 1592 #Meteor Histogram
1552 1593 finalHeights = finalMeteor[:,3]
@@ -1592,10 +1633,10 class WindProfiler(Operation):
1592 1633 def run(self, dataOut, technique, **kwargs):
1593 1634
1594 1635 param = dataOut.data_param
1595 if dataOut.abscissaList != None:
1596 absc = dataOut.abscissaList[:-1]
1636 # if dataOut.abscissaList != None:
1637 # absc = dataOut.abscissaList[:-1]
1597 1638 noise = dataOut.noise
1598 heightList = dataOut.getHeiRange()
1639 heightList = dataOut.heightList
1599 1640 SNR = dataOut.data_SNR
1600 1641
1601 1642 if technique == 'DBS':
@@ -1616,19 +1657,14 class WindProfiler(Operation):
1616 1657 else: correctFactor = 1
1617 1658 if kwargs.has_key('channelList'):
1618 1659 channelList = kwargs['channelList']
1660 if len(channelList) == 2:
1661 horizontalOnly = True
1619 1662 arrayChannel = numpy.array(channelList)
1620 1663 param = param[arrayChannel,:,:]
1621 1664 theta_x = theta_x[arrayChannel]
1622 1665 theta_y = theta_y[arrayChannel]
1623 1666
1624 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 1668 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
1633 1669 dataOut.utctimeInit = dataOut.utctime
1634 1670 dataOut.outputInterval = dataOut.timeInterval
@@ -1685,7 +1721,7 class WindProfiler(Operation):
1685 1721 if self.__isConfig == False:
1686 1722 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1687 1723 #Get Initial LTC time
1688 self.__initime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
1724 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1689 1725 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1690 1726
1691 1727 self.__isConfig = True
@@ -1695,7 +1731,7 class WindProfiler(Operation):
1695 1731 self.__firstdata = copy.copy(dataOut)
1696 1732
1697 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 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 1815 dataOut.outputInterval = dataOut.timeInterval
1780 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 2144 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now