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