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