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