##// END OF EJS Templates
schain
RSS
-Modifications to DBS module to accept azimuth and elevation angles as inputs....
Julio Valdez -
r503:ea980322a40d
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,579 +1,580
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from figure import Figure, isRealtime
6 6
7 7 class MomentsPlot(Figure):
8 8
9 9 isConfig = None
10 10 __nsubplots = None
11 11
12 12 WIDTHPROF = None
13 13 HEIGHTPROF = None
14 14 PREFIX = 'prm'
15 15
16 16 def __init__(self):
17 17
18 18 self.isConfig = False
19 19 self.__nsubplots = 1
20 20
21 21 self.WIDTH = 280
22 22 self.HEIGHT = 250
23 23 self.WIDTHPROF = 120
24 24 self.HEIGHTPROF = 0
25 25 self.counter_imagwr = 0
26 26
27 27 self.PLOT_CODE = 1
28 28 self.FTP_WEI = None
29 29 self.EXP_CODE = None
30 30 self.SUB_EXP_CODE = None
31 31 self.PLOT_POS = None
32 32
33 33 def getSubplots(self):
34 34
35 35 ncol = int(numpy.sqrt(self.nplots)+0.9)
36 36 nrow = int(self.nplots*1./ncol + 0.9)
37 37
38 38 return nrow, ncol
39 39
40 40 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
41 41
42 42 self.__showprofile = showprofile
43 43 self.nplots = nplots
44 44
45 45 ncolspan = 1
46 46 colspan = 1
47 47 if showprofile:
48 48 ncolspan = 3
49 49 colspan = 2
50 50 self.__nsubplots = 2
51 51
52 52 self.createFigure(id = id,
53 53 wintitle = wintitle,
54 54 widthplot = self.WIDTH + self.WIDTHPROF,
55 55 heightplot = self.HEIGHT + self.HEIGHTPROF,
56 56 show=show)
57 57
58 58 nrow, ncol = self.getSubplots()
59 59
60 60 counter = 0
61 61 for y in range(nrow):
62 62 for x in range(ncol):
63 63
64 64 if counter >= self.nplots:
65 65 break
66 66
67 67 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
68 68
69 69 if showprofile:
70 70 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
71 71
72 72 counter += 1
73 73
74 74 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
75 75 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
76 76 save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
77 77 server=None, folder=None, username=None, password=None,
78 78 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
79 79
80 80 """
81 81
82 82 Input:
83 83 dataOut :
84 84 id :
85 85 wintitle :
86 86 channelList :
87 87 showProfile :
88 88 xmin : None,
89 89 xmax : None,
90 90 ymin : None,
91 91 ymax : None,
92 92 zmin : None,
93 93 zmax : None
94 94 """
95 95
96 96 if dataOut.flagNoData:
97 97 return None
98 98
99 99 if realtime:
100 100 if not(isRealtime(utcdatatime = dataOut.utctime)):
101 101 print 'Skipping this plot function'
102 102 return
103 103
104 104 if channelList == None:
105 105 channelIndexList = dataOut.channelIndexList
106 106 else:
107 107 channelIndexList = []
108 108 for channel in channelList:
109 109 if channel not in dataOut.channelList:
110 110 raise ValueError, "Channel %d is not in dataOut.channelList"
111 111 channelIndexList.append(dataOut.channelList.index(channel))
112 112
113 113 factor = dataOut.normFactor
114 114 x = dataOut.abscissaRange
115 115 y = dataOut.heightRange
116 116
117 117 z = dataOut.data_pre[channelIndexList,:,:]/factor
118 118 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
119 119 avg = numpy.average(z, axis=1)
120 120 noise = dataOut.noise/factor
121 121
122 122 zdB = 10*numpy.log10(z)
123 123 avgdB = 10*numpy.log10(avg)
124 124 noisedB = 10*numpy.log10(noise)
125 125
126 126 #thisDatetime = dataOut.datatime
127 127 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
128 128 title = wintitle + " Parameters"
129 129 xlabel = "Velocity (m/s)"
130 130 ylabel = "Range (Km)"
131 131
132 132 if not self.isConfig:
133 133
134 134 nplots = len(channelIndexList)
135 135
136 136 self.setup(id=id,
137 137 nplots=nplots,
138 138 wintitle=wintitle,
139 139 showprofile=showprofile,
140 140 show=show)
141 141
142 142 if xmin == None: xmin = numpy.nanmin(x)
143 143 if xmax == None: xmax = numpy.nanmax(x)
144 144 if ymin == None: ymin = numpy.nanmin(y)
145 145 if ymax == None: ymax = numpy.nanmax(y)
146 146 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
147 147 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
148 148
149 149 self.FTP_WEI = ftp_wei
150 150 self.EXP_CODE = exp_code
151 151 self.SUB_EXP_CODE = sub_exp_code
152 152 self.PLOT_POS = plot_pos
153 153
154 154 self.isConfig = True
155 155
156 156 self.setWinTitle(title)
157 157
158 158 for i in range(self.nplots):
159 159 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
160 160 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
161 161 axes = self.axesList[i*self.__nsubplots]
162 162 axes.pcolor(x, y, zdB[i,:,:],
163 163 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
164 164 xlabel=xlabel, ylabel=ylabel, title=title,
165 165 ticksize=9, cblabel='')
166 166 #Mean Line
167 167 mean = dataOut.data_param[i, 1, :]
168 168 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
169 169
170 170 if self.__showprofile:
171 171 axes = self.axesList[i*self.__nsubplots +1]
172 172 axes.pline(avgdB[i], y,
173 173 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
174 174 xlabel='dB', ylabel='', title='',
175 175 ytick_visible=False,
176 176 grid='x')
177 177
178 178 noiseline = numpy.repeat(noisedB[i], len(y))
179 179 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
180 180
181 181 self.draw()
182 182
183 183 if figfile == None:
184 184 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
185 185 figfile = self.getFilename(name = str_datetime)
186 186
187 187 if figpath != '':
188 188 self.counter_imagwr += 1
189 189 if (self.counter_imagwr>=wr_period):
190 190 # store png plot to local folder
191 191 self.saveFigure(figpath, figfile)
192 192 # store png plot to FTP server according to RT-Web format
193 193 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
194 194 ftp_filename = os.path.join(figpath, name)
195 195 self.saveFigure(figpath, ftp_filename)
196 196 self.counter_imagwr = 0
197 197
198 198 class SkyMapPlot(Figure):
199 199
200 200 __isConfig = None
201 201 __nsubplots = None
202 202
203 203 WIDTHPROF = None
204 204 HEIGHTPROF = None
205 205 PREFIX = 'prm'
206 206
207 207 def __init__(self):
208 208
209 209 self.__isConfig = False
210 210 self.__nsubplots = 1
211 211
212 212 # self.WIDTH = 280
213 213 # self.HEIGHT = 250
214 214 self.WIDTH = 600
215 215 self.HEIGHT = 600
216 216 self.WIDTHPROF = 120
217 217 self.HEIGHTPROF = 0
218 218 self.counter_imagwr = 0
219 219
220 220 self.PLOT_CODE = 1
221 221 self.FTP_WEI = None
222 222 self.EXP_CODE = None
223 223 self.SUB_EXP_CODE = None
224 224 self.PLOT_POS = None
225 225
226 226 def getSubplots(self):
227 227
228 228 ncol = int(numpy.sqrt(self.nplots)+0.9)
229 229 nrow = int(self.nplots*1./ncol + 0.9)
230 230
231 231 return nrow, ncol
232 232
233 233 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
234 234
235 235 self.__showprofile = showprofile
236 236 self.nplots = nplots
237 237
238 238 ncolspan = 1
239 239 colspan = 1
240 240
241 241 self.createFigure(id = id,
242 242 wintitle = wintitle,
243 243 widthplot = self.WIDTH, #+ self.WIDTHPROF,
244 244 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
245 245 show=show)
246 246
247 247 nrow, ncol = 1,1
248 248 counter = 0
249 249 x = 0
250 250 y = 0
251 251 self.addAxes(1, 1, 0, 0, 1, 1, True)
252 252
253 253 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
254 254 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
255 255 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
256 256 server=None, folder=None, username=None, password=None,
257 257 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
258 258
259 259 """
260 260
261 261 Input:
262 262 dataOut :
263 263 id :
264 264 wintitle :
265 265 channelList :
266 266 showProfile :
267 267 xmin : None,
268 268 xmax : None,
269 269 ymin : None,
270 270 ymax : None,
271 271 zmin : None,
272 272 zmax : None
273 273 """
274 274
275 275 arrayParameters = dataOut.data_param
276 276 error = arrayParameters[:,-1]
277 277 indValid = numpy.where(error == 0)[0]
278 278 finalMeteor = arrayParameters[indValid,:]
279 279 finalAzimuth = finalMeteor[:,4]
280 280 finalZenith = finalMeteor[:,5]
281 281
282 282 x = finalAzimuth*numpy.pi/180
283 283 y = finalZenith
284 284
285 285
286 286 #thisDatetime = dataOut.datatime
287 287 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
288 288 title = wintitle + " Parameters"
289 289 xlabel = "Zonal Zenith Angle (deg) "
290 290 ylabel = "Meridional Zenith Angle (deg)"
291 291
292 292 if not self.__isConfig:
293 293
294 294 nplots = 1
295 295
296 296 self.setup(id=id,
297 297 nplots=nplots,
298 298 wintitle=wintitle,
299 299 showprofile=showprofile,
300 300 show=show)
301 301
302 302 self.FTP_WEI = ftp_wei
303 303 self.EXP_CODE = exp_code
304 304 self.SUB_EXP_CODE = sub_exp_code
305 305 self.PLOT_POS = plot_pos
306 306 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
307 307 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
308 308 self.__isConfig = True
309 309
310 310 self.setWinTitle(title)
311 311
312 312 i = 0
313 313 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
314 314
315 315 axes = self.axesList[i*self.__nsubplots]
316 316 nevents = axes.x_buffer.shape[0] + x.shape[0]
317 317 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
318 318 axes.polar(x, y,
319 319 title=title, xlabel=xlabel, ylabel=ylabel,
320 320 ticksize=9, cblabel='')
321 321
322 322 self.draw()
323 323
324 324 if save:
325 325
326 326 self.counter_imagwr += 1
327 327 if (self.counter_imagwr==wr_period):
328 328
329 329 if figfile == None:
330 330 figfile = self.getFilename(name = self.name)
331 331 self.saveFigure(figpath, figfile)
332 332
333 333 if ftp:
334 334 #provisionalmente envia archivos en el formato de la web en tiempo real
335 335 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
336 336 path = '%s%03d' %(self.PREFIX, self.id)
337 337 ftp_file = os.path.join(path,'ftp','%s.png'%name)
338 338 self.saveFigure(figpath, ftp_file)
339 339 ftp_filename = os.path.join(figpath,ftp_file)
340 340
341 341
342 342 try:
343 343 self.sendByFTP(ftp_filename, server, folder, username, password)
344 344 except:
345 345 self.counter_imagwr = 0
346 346 raise ValueError, 'Error FTP'
347 347
348 348 self.counter_imagwr = 0
349 349
350 350
351 351 class WindProfilerPlot(Figure):
352 352
353 353 __isConfig = None
354 354 __nsubplots = None
355 355
356 356 WIDTHPROF = None
357 357 HEIGHTPROF = None
358 358 PREFIX = 'wind'
359 359
360 360 def __init__(self):
361 361
362 362 self.timerange = 2*60*60
363 363 self.isConfig = False
364 364 self.__nsubplots = 1
365 365
366 366 self.WIDTH = 800
367 367 self.HEIGHT = 150
368 368 self.WIDTHPROF = 120
369 369 self.HEIGHTPROF = 0
370 370 self.counter_imagwr = 0
371 371
372 372 self.PLOT_CODE = 0
373 373 self.FTP_WEI = None
374 374 self.EXP_CODE = None
375 375 self.SUB_EXP_CODE = None
376 376 self.PLOT_POS = None
377 377 self.tmin = None
378 378 self.tmax = None
379 379
380 380 self.xmin = None
381 381 self.xmax = None
382 382
383 383 self.figfile = None
384 384
385 385 def getSubplots(self):
386 386
387 387 ncol = 1
388 388 nrow = self.nplots
389 389
390 390 return nrow, ncol
391 391
392 392 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
393 393
394 394 self.__showprofile = showprofile
395 395 self.nplots = nplots
396 396
397 397 ncolspan = 1
398 398 colspan = 1
399 399
400 400 self.createFigure(id = id,
401 401 wintitle = wintitle,
402 402 widthplot = self.WIDTH + self.WIDTHPROF,
403 403 heightplot = self.HEIGHT + self.HEIGHTPROF,
404 404 show=show)
405 405
406 406 nrow, ncol = self.getSubplots()
407 407
408 408 counter = 0
409 409 for y in range(nrow):
410 410 if counter >= self.nplots:
411 411 break
412 412
413 413 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
414 414 counter += 1
415 415
416 416 def run(self, dataOut, id, wintitle="", channelList=None,
417 417 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
418 418 zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
419 419 timerange=None, SNRthresh = None,
420 420 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
421 421 server=None, folder=None, username=None, password=None,
422 422 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
423 423 """
424 424
425 425 Input:
426 426 dataOut :
427 427 id :
428 428 wintitle :
429 429 channelList :
430 430 showProfile :
431 431 xmin : None,
432 432 xmax : None,
433 433 ymin : None,
434 434 ymax : None,
435 435 zmin : None,
436 436 zmax : None
437 437 """
438 438
439 439 if channelList == None:
440 440 channelIndexList = dataOut.channelIndexList
441 441 else:
442 442 channelIndexList = []
443 443 for channel in channelList:
444 444 if channel not in dataOut.channelList:
445 445 raise ValueError, "Channel %d is not in dataOut.channelList"
446 446 channelIndexList.append(dataOut.channelList.index(channel))
447 447
448 448 if timerange != None:
449 449 self.timerange = timerange
450 450
451 451 tmin = None
452 452 tmax = None
453 453
454 454 x = dataOut.getTimeRange1()
455 # y = dataOut.heightRange
455 456 y = dataOut.heightRange
456 457
457 458 z = dataOut.winds
458 459 nplots = z.shape[0] #Number of wind dimensions estimated
459 460 nplotsw = nplots
460 461
461 462 #If there is a SNR function defined
462 463 if dataOut.SNR != None:
463 464 nplots += 1
464 465 SNR = dataOut.SNR
465 466 SNRavg = numpy.average(SNR, axis=0)
466 467
467 468 SNRdB = 10*numpy.log10(SNR)
468 469 SNRavgdB = 10*numpy.log10(SNRavg)
469 470
470 471 if SNRthresh == None: SNRthresh = -5.0
471 472 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
472 473
473 474 for i in range(nplotsw):
474 475 z[i,ind] = numpy.nan
475 476
476 477
477 478 showprofile = False
478 479 # thisDatetime = dataOut.datatime
479 480 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
480 481 title = wintitle + "Wind"
481 482 xlabel = ""
482 483 ylabel = "Range (Km)"
483 484
484 485 if not self.__isConfig:
485 486
486 487
487 488
488 489 self.setup(id=id,
489 490 nplots=nplots,
490 491 wintitle=wintitle,
491 492 showprofile=showprofile,
492 493 show=show)
493 494
494 495 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
495 496
496 497 if ymin == None: ymin = numpy.nanmin(y)
497 498 if ymax == None: ymax = numpy.nanmax(y)
498 499
499 500 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
500 501 #if numpy.isnan(zmax): zmax = 50
501 502 if zmin == None: zmin = -zmax
502 503
503 504 if nplotsw == 3:
504 505 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
505 506 if zmin_ver == None: zmin_ver = -zmax_ver
506 507
507 508 if dataOut.SNR != None:
508 509 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
509 510 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
510 511
511 512 self.FTP_WEI = ftp_wei
512 513 self.EXP_CODE = exp_code
513 514 self.SUB_EXP_CODE = sub_exp_code
514 515 self.PLOT_POS = plot_pos
515 516
516 517 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
517 518 self.__isConfig = True
518 519
519 520
520 521 self.setWinTitle(title)
521 522
522 523 if ((self.xmax - x[1]) < (x[1]-x[0])):
523 524 x[1] = self.xmax
524 525
525 526 strWind = ['Zonal', 'Meridional', 'Vertical']
526 527 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
527 528 zmaxVector = [zmax, zmax, zmax_ver]
528 529 zminVector = [zmin, zmin, zmin_ver]
529 530 windFactor = [1,1,100]
530 531
531 532 for i in range(nplotsw):
532 533
533 534 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
534 535 axes = self.axesList[i*self.__nsubplots]
535 536
536 537 z1 = z[i,:].reshape((1,-1))*windFactor[i]
537 538
538 539 axes.pcolorbuffer(x, y, z1,
539 540 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
540 541 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
541 542 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
542 543
543 544 if dataOut.SNR != None:
544 545 i += 1
545 546 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
546 547 axes = self.axesList[i*self.__nsubplots]
547 548
548 549 SNRavgdB = SNRavgdB.reshape((1,-1))
549 550
550 551 axes.pcolorbuffer(x, y, SNRavgdB,
551 552 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
552 553 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
553 554 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
554 555
555 556 self.draw()
556 557
557 558 if x[1] >= self.axesList[0].xmax:
558 559 self.counter_imagwr = wr_period
559 560 self.__isConfig = False
560 561
561 562
562 563 if self.figfile == None:
563 564 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
564 565 self.figfile = self.getFilename(name = str_datetime)
565 566
566 567 if figpath != '':
567 568
568 569 self.counter_imagwr += 1
569 570 if (self.counter_imagwr>=wr_period):
570 571 # store png plot to local folder
571 572 self.saveFigure(figpath, self.figfile)
572 573 # store png plot to FTP server according to RT-Web format
573 574 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
574 575 ftp_filename = os.path.join(figpath, name)
575 576 self.saveFigure(figpath, ftp_filename)
576 577
577 578 self.counter_imagwr = 0
578 579
579 580 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1521 +1,1539
1 1 import numpy
2 2 import math
3 3 from scipy import optimize
4 4 from scipy import interpolate
5 5 from scipy import signal
6 6 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10 10
11 11
12 12 from jroproc_base import ProcessingUnit, Operation
13 13 from model.data.jrodata import Parameters
14 14
15 15
16 16 class ParametersProc(ProcessingUnit):
17 17
18 18 nSeconds = None
19 19
20 20 def __init__(self):
21 21 ProcessingUnit.__init__(self)
22 22
23 23 self.objectDict = {}
24 24 self.buffer = None
25 25 self.firstdatatime = None
26 26 self.profIndex = 0
27 27 self.dataOut = Parameters()
28 28
29 29 def __updateObjFromInput(self):
30 30
31 31 self.dataOut.inputUnit = self.dataIn.type
32 32
33 33 self.dataOut.timeZone = self.dataIn.timeZone
34 34 self.dataOut.dstFlag = self.dataIn.dstFlag
35 35 self.dataOut.errorCount = self.dataIn.errorCount
36 36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37 37
38 38 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
39 39 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
40 40 self.dataOut.channelList = self.dataIn.channelList
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
43 43 # self.dataOut.nHeights = self.dataIn.nHeights
44 44 # self.dataOut.nChannels = self.dataIn.nChannels
45 45 self.dataOut.nBaud = self.dataIn.nBaud
46 46 self.dataOut.nCode = self.dataIn.nCode
47 47 self.dataOut.code = self.dataIn.code
48 48 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
49 49 self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
50 50 self.dataOut.utctime = self.firstdatatime
51 51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 53 # self.dataOut.nCohInt = self.dataIn.nCohInt
54 54 # self.dataOut.nIncohInt = 1
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 58 self.dataOut.heightRange = self.dataIn.getHeiRange()
59 59 self.dataOut.frequency = self.dataIn.frequency
60 60
61 61 def run(self, nSeconds = None, nProfiles = None):
62 62
63 63 self.dataOut.flagNoData = True
64 64
65 65 if self.firstdatatime == None:
66 66 self.firstdatatime = self.dataIn.utctime
67 67
68 68 #---------------------- Voltage Data ---------------------------
69 69
70 70 if self.dataIn.type == "Voltage":
71 71 if nSeconds != None:
72 72 self.nSeconds = nSeconds
73 73 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
74 74
75 75 if self.buffer == None:
76 76 self.buffer = numpy.zeros((self.dataIn.nChannels,
77 77 self.nProfiles,
78 78 self.dataIn.nHeights),
79 79 dtype='complex')
80 80
81 81 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
82 82 self.profIndex += 1
83 83
84 84 if self.profIndex == self.nProfiles:
85 85
86 86 self.__updateObjFromInput()
87 87 self.dataOut.data_pre = self.buffer.copy()
88 88 self.dataOut.paramInterval = nSeconds
89 89 self.dataOut.flagNoData = False
90 90
91 91 self.buffer = None
92 92 self.firstdatatime = None
93 93 self.profIndex = 0
94 94
95 95 #---------------------- Spectra Data ---------------------------
96 96
97 97 if self.dataIn.type == "Spectra":
98 98 self.dataOut.data_pre = self.dataIn.data_spc.copy()
99 99 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
100 100 self.dataOut.noise = self.dataIn.getNoise()
101 101 self.dataOut.normFactor = self.dataIn.normFactor
102 102
103 103 self.__updateObjFromInput()
104 104 self.dataOut.flagNoData = False
105 105 self.firstdatatime = None
106 106
107 107 #---------------------- Correlation Data ---------------------------
108 108
109 109 if self.dataIn.type == "Correlation":
110 110 lagRRange = self.dataIn.lagR
111 111 indR = numpy.where(lagRRange == 0)[0][0]
112 112
113 113 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
114 114 self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
115 115 self.dataOut.noise = self.dataIn.noise
116 116 self.dataOut.normFactor = self.dataIn.normFactor
117 117 self.dataOut.SNR = self.dataIn.SNR
118 118 self.dataOut.pairsList = self.dataIn.pairsList
119 119
120 120 self.__updateObjFromInput()
121 121 self.dataOut.flagNoData = False
122 122 self.firstdatatime = None
123 123
124 124 #------------------- Get Moments ----------------------------------
125 125 def GetMoments(self, channelList = None):
126 126 '''
127 127 Function GetMoments()
128 128
129 129 Input:
130 130 channelList : simple channel list to select e.g. [2,3,7]
131 131 self.dataOut.data_pre
132 132 self.dataOut.abscissaRange
133 133 self.dataOut.noise
134 134
135 135 Affected:
136 136 self.dataOut.data_param
137 137 self.dataOut.SNR
138 138
139 139 '''
140 140 data = self.dataOut.data_pre
141 141 absc = self.dataOut.abscissaRange[:-1]
142 142 noise = self.dataOut.noise
143 143
144 144 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
145 145
146 146 if channelList== None: channelList = self.dataOut.channelList
147 147
148 148 for ind in channelList:
149 149 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
150 150
151 151 self.dataOut.data_param = data_param[:,1:]
152 152 self.dataOut.SNR = data_param[:,0]
153 153 return
154 154
155 155 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):
156 156
157 157 if (nicoh == None): nicoh = 1
158 158 if (graph == None): graph = 0
159 159 if (smooth == None): smooth = 0
160 160 elif (self.smooth < 3): smooth = 0
161 161
162 162 if (type1 == None): type1 = 0
163 163 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
164 164 if (snrth == None): snrth = -3
165 165 if (dc == None): dc = 0
166 166 if (aliasing == None): aliasing = 0
167 167 if (oldfd == None): oldfd = 0
168 168 if (wwauto == None): wwauto = 0
169 169
170 170 if (n0 < 1.e-20): n0 = 1.e-20
171 171
172 172 freq = oldfreq
173 173 vec_power = numpy.zeros(oldspec.shape[1])
174 174 vec_fd = numpy.zeros(oldspec.shape[1])
175 175 vec_w = numpy.zeros(oldspec.shape[1])
176 176 vec_snr = numpy.zeros(oldspec.shape[1])
177 177
178 178 for ind in range(oldspec.shape[1]):
179 179
180 180 spec = oldspec[:,ind]
181 181 aux = spec*fwindow
182 182 max_spec = aux.max()
183 183 m = list(aux).index(max_spec)
184 184
185 185 #Smooth
186 186 if (smooth == 0): spec2 = spec
187 187 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
188 188
189 189 # Calculo de Momentos
190 190 bb = spec2[range(m,spec2.size)]
191 191 bb = (bb<n0).nonzero()
192 192 bb = bb[0]
193 193
194 194 ss = spec2[range(0,m + 1)]
195 195 ss = (ss<n0).nonzero()
196 196 ss = ss[0]
197 197
198 198 if (bb.size == 0):
199 199 bb0 = spec.size - 1 - m
200 200 else:
201 201 bb0 = bb[0] - 1
202 202 if (bb0 < 0):
203 203 bb0 = 0
204 204
205 205 if (ss.size == 0): ss1 = 1
206 206 else: ss1 = max(ss) + 1
207 207
208 208 if (ss1 > m): ss1 = m
209 209
210 210 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
211 211 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
212 212 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
213 213 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
214 214 snr = (spec2.mean()-n0)/n0
215 215
216 216 if (snr < 1.e-20) :
217 217 snr = 1.e-20
218 218
219 219 vec_power[ind] = power
220 220 vec_fd[ind] = fd
221 221 vec_w[ind] = w
222 222 vec_snr[ind] = snr
223 223
224 224 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
225 225 return moments
226 226
227 227 #------------------- Get Lags ----------------------------------
228 228
229 229 def GetLags(self):
230 230 '''
231 231 Function GetMoments()
232 232
233 233 Input:
234 234 self.dataOut.data_pre
235 235 self.dataOut.abscissaRange
236 236 self.dataOut.noise
237 237 self.dataOut.normFactor
238 238 self.dataOut.SNR
239 239 self.dataOut.pairsList
240 240 self.dataOut.nChannels
241 241
242 242 Affected:
243 243 self.dataOut.data_param
244 244
245 245 '''
246 246 data = self.dataOut.data_pre
247 247 normFactor = self.dataOut.normFactor
248 248 nHeights = self.dataOut.nHeights
249 249 absc = self.dataOut.abscissaRange[:-1]
250 250 noise = self.dataOut.noise
251 251 SNR = self.dataOut.SNR
252 252 pairsList = self.dataOut.pairsList
253 253 nChannels = self.dataOut.nChannels
254 254 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
255 255 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
256 256
257 257 dataNorm = numpy.abs(data)
258 258 for l in range(len(pairsList)):
259 259 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
260 260
261 261 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
262 262 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
263 263 return
264 264
265 265 def __getPairsAutoCorr(self, pairsList, nChannels):
266 266
267 267 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
268 268
269 269 for l in range(len(pairsList)):
270 270 firstChannel = pairsList[l][0]
271 271 secondChannel = pairsList[l][1]
272 272
273 273 #Obteniendo pares de Autocorrelacion
274 274 if firstChannel == secondChannel:
275 275 pairsAutoCorr[firstChannel] = int(l)
276 276
277 277 pairsAutoCorr = pairsAutoCorr.astype(int)
278 278
279 279 pairsCrossCorr = range(len(pairsList))
280 280 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
281 281
282 282 return pairsAutoCorr, pairsCrossCorr
283 283
284 284 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
285 285
286 286 Pt0 = data.shape[1]/2
287 287 #Funcion de Autocorrelacion
288 288 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
289 289
290 290 #Obtencion Indice de TauCross
291 291 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
292 292 #Obtencion Indice de TauAuto
293 293 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
294 294 CCValue = data[pairsCrossCorr,Pt0,:]
295 295 for i in range(pairsCrossCorr.size):
296 296 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
297 297
298 298 #Obtencion de TauCross y TauAuto
299 299 tauCross = lagTRange[indCross]
300 300 tauAuto = lagTRange[indAuto]
301 301
302 302 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
303 303
304 304 tauCross[Nan1,Nan2] = numpy.nan
305 305 tauAuto[Nan1,Nan2] = numpy.nan
306 306 tau = numpy.vstack((tauCross,tauAuto))
307 307
308 308 return tau
309 309
310 310 def __calculateLag1Phase(self, data, pairs, lagTRange):
311 311 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
312 312 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
313 313
314 314 phase = numpy.angle(data1[lag1,:])
315 315
316 316 return phase
317 317 #------------------- Detect Meteors ------------------------------
318 318
319 319 def DetectMeteors(self, hei_ref = None, tauindex = 0,
320 320 predefinedPhaseShifts = None, centerReceiverIndex = 2,
321 321 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
322 322 noise_timeStep = 4, noise_multiple = 4,
323 323 multDet_timeLimit = 1, multDet_rangeLimit = 3,
324 324 phaseThresh = 20, SNRThresh = 8,
325 325 hmin = 70, hmax=110, azimuth = 0) :
326 326
327 327 '''
328 328 Function DetectMeteors()
329 329 Project developed with paper:
330 330 HOLDSWORTH ET AL. 2004
331 331
332 332 Input:
333 333 self.dataOut.data_pre
334 334
335 335 centerReceiverIndex: From the channels, which is the center receiver
336 336
337 337 hei_ref: Height reference for the Beacon signal extraction
338 338 tauindex:
339 339 predefinedPhaseShifts: Predefined phase offset for the voltge signals
340 340
341 341 cohDetection: Whether to user Coherent detection or not
342 342 cohDet_timeStep: Coherent Detection calculation time step
343 343 cohDet_thresh: Coherent Detection phase threshold to correct phases
344 344
345 345 noise_timeStep: Noise calculation time step
346 346 noise_multiple: Noise multiple to define signal threshold
347 347
348 348 multDet_timeLimit: Multiple Detection Removal time limit in seconds
349 349 multDet_rangeLimit: Multiple Detection Removal range limit in km
350 350
351 351 phaseThresh: Maximum phase difference between receiver to be consider a meteor
352 352 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
353 353
354 354 hmin: Minimum Height of the meteor to use it in the further wind estimations
355 355 hmax: Maximum Height of the meteor to use it in the further wind estimations
356 356 azimuth: Azimuth angle correction
357 357
358 358 Affected:
359 359 self.dataOut.data_param
360 360
361 361 Rejection Criteria (Errors):
362 362 0: No error; analysis OK
363 363 1: SNR < SNR threshold
364 364 2: angle of arrival (AOA) ambiguously determined
365 365 3: AOA estimate not feasible
366 366 4: Large difference in AOAs obtained from different antenna baselines
367 367 5: echo at start or end of time series
368 368 6: echo less than 5 examples long; too short for analysis
369 369 7: echo rise exceeds 0.3s
370 370 8: echo decay time less than twice rise time
371 371 9: large power level before echo
372 372 10: large power level after echo
373 373 11: poor fit to amplitude for estimation of decay time
374 374 12: poor fit to CCF phase variation for estimation of radial drift velocity
375 375 13: height unresolvable echo: not valid height within 70 to 110 km
376 376 14: height ambiguous echo: more then one possible height within 70 to 110 km
377 377 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
378 378 16: oscilatory echo, indicating event most likely not an underdense echo
379 379
380 380 17: phase difference in meteor Reestimation
381 381
382 382 Data Storage:
383 383 Meteors for Wind Estimation (8):
384 384 Day Hour | Range Height
385 385 Azimuth Zenith errorCosDir
386 386 VelRad errorVelRad
387 387 TypeError
388 388
389 389 '''
390 390 #Get Beacon signal
391 391 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
392 392
393 393 if hei_ref != None:
394 394 newheis = numpy.where(self.dataOut.heightList>hei_ref)
395 395
396 396 heiRang = self.dataOut.getHeiRange()
397 397 #Pairs List
398 398 pairslist = []
399 399 nChannel = self.dataOut.nChannels
400 400 for i in range(nChannel):
401 401 if i != centerReceiverIndex:
402 402 pairslist.append((centerReceiverIndex,i))
403 403
404 404 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
405 405 # see if the user put in pre defined phase shifts
406 406 voltsPShift = self.dataOut.data_pre.copy()
407 407
408 408 if predefinedPhaseShifts != None:
409 409 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
410 410 else:
411 411 #get hardware phase shifts using beacon signal
412 412 hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
413 413 hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
414 414
415 415 voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
416 416 for i in range(self.dataOut.data_pre.shape[0]):
417 417 voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
418 418 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
419 419
420 420 #Remove DC
421 421 voltsDC = numpy.mean(voltsPShift,1)
422 422 voltsDC = numpy.mean(voltsDC,1)
423 423 for i in range(voltsDC.shape[0]):
424 424 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
425 425
426 426 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
427 427 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
428 428
429 429 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
430 430 #Coherent Detection
431 431 if cohDetection:
432 432 #use coherent detection to get the net power
433 433 cohDet_thresh = cohDet_thresh*numpy.pi/180
434 434 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
435 435
436 436 #Non-coherent detection!
437 437 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
438 438 #********** END OF COH/NON-COH POWER CALCULATION**********************
439 439
440 440 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
441 441 #Get noise
442 442 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
443 443 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
444 444 #Get signal threshold
445 445 signalThresh = noise_multiple*noise
446 446 #Meteor echoes detection
447 447 listMeteors = self.__findMeteors(powerNet, signalThresh)
448 448 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
449 449
450 450 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
451 451 #Parameters
452 452 heiRange = self.dataOut.getHeiRange()
453 453 rangeInterval = heiRange[1] - heiRange[0]
454 454 rangeLimit = multDet_rangeLimit/rangeInterval
455 455 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
456 456 #Multiple detection removals
457 457 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
458 458 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
459 459
460 460 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
461 461 #Parameters
462 462 phaseThresh = phaseThresh*numpy.pi/180
463 463 thresh = [phaseThresh, noise_multiple, SNRThresh]
464 464 #Meteor reestimation (Errors N 1, 6, 12, 17)
465 465 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
466 466 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
467 467 #Estimation of decay times (Errors N 7, 8, 11)
468 468 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
469 469 #******************* END OF METEOR REESTIMATION *******************
470 470
471 471 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
472 472 #Calculating Radial Velocity (Error N 15)
473 473 radialStdThresh = 10
474 474 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
475 475
476 476 if len(listMeteors4) > 0:
477 477 #Setting New Array
478 478 date = repr(self.dataOut.datatime)
479 479 arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
480 480
481 481 #Calculate AOA (Error N 3, 4)
482 482 #JONES ET AL. 1998
483 483 AOAthresh = numpy.pi/8
484 484 error = arrayParameters[:,-1]
485 485 phases = -arrayMeteors4[:,9:13]
486 486 pairsList = []
487 487 pairsList.append((0,3))
488 488 pairsList.append((1,2))
489 489 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
490 490
491 491 #Calculate Heights (Error N 13 and 14)
492 492 error = arrayParameters[:,-1]
493 493 Ranges = arrayParameters[:,2]
494 494 zenith = arrayParameters[:,5]
495 495 arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
496 496 #********************* END OF PARAMETERS CALCULATION **************************
497 497
498 498 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
499 499 self.dataOut.data_param = arrayParameters
500 500
501 501 return
502 502
503 503 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
504 504
505 505 minIndex = min(newheis[0])
506 506 maxIndex = max(newheis[0])
507 507
508 508 voltage = voltage0[:,:,minIndex:maxIndex+1]
509 509 nLength = voltage.shape[1]/n
510 510 nMin = 0
511 511 nMax = 0
512 512 phaseOffset = numpy.zeros((len(pairslist),n))
513 513
514 514 for i in range(n):
515 515 nMax += nLength
516 516 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
517 517 phaseCCF = numpy.mean(phaseCCF, axis = 2)
518 518 phaseOffset[:,i] = phaseCCF.transpose()
519 519 nMin = nMax
520 520 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
521 521
522 522 #Remove Outliers
523 523 factor = 2
524 524 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
525 525 dw = numpy.std(wt,axis = 1)
526 526 dw = dw.reshape((dw.size,1))
527 527 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
528 528 phaseOffset[ind] = numpy.nan
529 529 phaseOffset = stats.nanmean(phaseOffset, axis=1)
530 530
531 531 return phaseOffset
532 532
533 533 def __shiftPhase(self, data, phaseShift):
534 534 #this will shift the phase of a complex number
535 535 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
536 536 return dataShifted
537 537
538 538 def __estimatePhaseDifference(self, array, pairslist):
539 539 nChannel = array.shape[0]
540 540 nHeights = array.shape[2]
541 541 numPairs = len(pairslist)
542 542 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
543 543 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
544 544
545 545 #Correct phases
546 546 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
547 547 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
548 548
549 549 if indDer[0].shape[0] > 0:
550 550 for i in range(indDer[0].shape[0]):
551 551 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
552 552 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
553 553
554 554 # for j in range(numSides):
555 555 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
556 556 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
557 557 #
558 558 #Linear
559 559 phaseInt = numpy.zeros((numPairs,1))
560 560 angAllCCF = phaseCCF[:,[0,1,3,4],0]
561 561 for j in range(numPairs):
562 562 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
563 563 phaseInt[j] = fit[1]
564 564 #Phase Differences
565 565 phaseDiff = phaseInt - phaseCCF[:,2,:]
566 566 phaseArrival = phaseInt.reshape(phaseInt.size)
567 567
568 568 #Dealias
569 569 indAlias = numpy.where(phaseArrival > numpy.pi)
570 570 phaseArrival[indAlias] -= 2*numpy.pi
571 571 indAlias = numpy.where(phaseArrival < -numpy.pi)
572 572 phaseArrival[indAlias] += 2*numpy.pi
573 573
574 574 return phaseDiff, phaseArrival
575 575
576 576 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
577 577 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
578 578 #find the phase shifts of each channel over 1 second intervals
579 579 #only look at ranges below the beacon signal
580 580 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
581 581 numBlocks = int(volts.shape[1]/numProfPerBlock)
582 582 numHeights = volts.shape[2]
583 583 nChannel = volts.shape[0]
584 584 voltsCohDet = volts.copy()
585 585
586 586 pairsarray = numpy.array(pairslist)
587 587 indSides = pairsarray[:,1]
588 588 # indSides = numpy.array(range(nChannel))
589 589 # indSides = numpy.delete(indSides, indCenter)
590 590 #
591 591 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
592 592 listBlocks = numpy.array_split(volts, numBlocks, 1)
593 593
594 594 startInd = 0
595 595 endInd = 0
596 596
597 597 for i in range(numBlocks):
598 598 startInd = endInd
599 599 endInd = endInd + listBlocks[i].shape[1]
600 600
601 601 arrayBlock = listBlocks[i]
602 602 # arrayBlockCenter = listCenter[i]
603 603
604 604 #Estimate the Phase Difference
605 605 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
606 606 #Phase Difference RMS
607 607 arrayPhaseRMS = numpy.abs(phaseDiff)
608 608 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
609 609 indPhase = numpy.where(phaseRMSaux==4)
610 610 #Shifting
611 611 if indPhase[0].shape[0] > 0:
612 612 for j in range(indSides.size):
613 613 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
614 614 voltsCohDet[:,startInd:endInd,:] = arrayBlock
615 615
616 616 return voltsCohDet
617 617
618 618 def __calculateCCF(self, volts, pairslist ,laglist):
619 619
620 620 nHeights = volts.shape[2]
621 621 nPoints = volts.shape[1]
622 622 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
623 623
624 624 for i in range(len(pairslist)):
625 625 volts1 = volts[pairslist[i][0]]
626 626 volts2 = volts[pairslist[i][1]]
627 627
628 628 for t in range(len(laglist)):
629 629 idxT = laglist[t]
630 630 if idxT >= 0:
631 631 vStacked = numpy.vstack((volts2[idxT:,:],
632 632 numpy.zeros((idxT, nHeights),dtype='complex')))
633 633 else:
634 634 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
635 635 volts2[:(nPoints + idxT),:]))
636 636 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
637 637
638 638 vStacked = None
639 639 return voltsCCF
640 640
641 641 def __getNoise(self, power, timeSegment, timeInterval):
642 642 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
643 643 numBlocks = int(power.shape[0]/numProfPerBlock)
644 644 numHeights = power.shape[1]
645 645
646 646 listPower = numpy.array_split(power, numBlocks, 0)
647 647 noise = numpy.zeros((power.shape[0], power.shape[1]))
648 648 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
649 649
650 650 startInd = 0
651 651 endInd = 0
652 652
653 653 for i in range(numBlocks): #split por canal
654 654 startInd = endInd
655 655 endInd = endInd + listPower[i].shape[0]
656 656
657 657 arrayBlock = listPower[i]
658 658 noiseAux = numpy.mean(arrayBlock, 0)
659 659 # noiseAux = numpy.median(noiseAux)
660 660 # noiseAux = numpy.mean(arrayBlock)
661 661 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
662 662
663 663 noiseAux1 = numpy.mean(arrayBlock)
664 664 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
665 665
666 666 return noise, noise1
667 667
668 668 def __findMeteors(self, power, thresh):
669 669 nProf = power.shape[0]
670 670 nHeights = power.shape[1]
671 671 listMeteors = []
672 672
673 673 for i in range(nHeights):
674 674 powerAux = power[:,i]
675 675 threshAux = thresh[:,i]
676 676
677 677 indUPthresh = numpy.where(powerAux > threshAux)[0]
678 678 indDNthresh = numpy.where(powerAux <= threshAux)[0]
679 679
680 680 j = 0
681 681
682 682 while (j < indUPthresh.size - 2):
683 683 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
684 684 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
685 685 indDNthresh = indDNthresh[indDNAux]
686 686
687 687 if (indDNthresh.size > 0):
688 688 indEnd = indDNthresh[0] - 1
689 689 indInit = indUPthresh[j]
690 690
691 691 meteor = powerAux[indInit:indEnd + 1]
692 692 indPeak = meteor.argmax() + indInit
693 693 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
694 694
695 695 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
696 696 j = numpy.where(indUPthresh == indEnd)[0] + 1
697 697 else: j+=1
698 698 else: j+=1
699 699
700 700 return listMeteors
701 701
702 702 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
703 703
704 704 arrayMeteors = numpy.asarray(listMeteors)
705 705 listMeteors1 = []
706 706
707 707 while arrayMeteors.shape[0] > 0:
708 708 FLAs = arrayMeteors[:,4]
709 709 maxFLA = FLAs.argmax()
710 710 listMeteors1.append(arrayMeteors[maxFLA,:])
711 711
712 712 MeteorInitTime = arrayMeteors[maxFLA,1]
713 713 MeteorEndTime = arrayMeteors[maxFLA,3]
714 714 MeteorHeight = arrayMeteors[maxFLA,0]
715 715
716 716 #Check neighborhood
717 717 maxHeightIndex = MeteorHeight + rangeLimit
718 718 minHeightIndex = MeteorHeight - rangeLimit
719 719 minTimeIndex = MeteorInitTime - timeLimit
720 720 maxTimeIndex = MeteorEndTime + timeLimit
721 721
722 722 #Check Heights
723 723 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
724 724 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
725 725 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
726 726
727 727 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
728 728
729 729 return listMeteors1
730 730
731 731 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
732 732 numHeights = volts.shape[2]
733 733 nChannel = volts.shape[0]
734 734
735 735 thresholdPhase = thresh[0]
736 736 thresholdNoise = thresh[1]
737 737 thresholdDB = float(thresh[2])
738 738
739 739 thresholdDB1 = 10**(thresholdDB/10)
740 740 pairsarray = numpy.array(pairslist)
741 741 indSides = pairsarray[:,1]
742 742
743 743 pairslist1 = list(pairslist)
744 744 pairslist1.append((0,1))
745 745 pairslist1.append((3,4))
746 746
747 747 listMeteors1 = []
748 748 listPowerSeries = []
749 749 listVoltageSeries = []
750 750 #volts has the war data
751 751
752 752 if frequency == 30e6:
753 753 timeLag = 45*10**-3
754 754 else:
755 755 timeLag = 15*10**-3
756 756 lag = numpy.ceil(timeLag/timeInterval)
757 757
758 758 for i in range(len(listMeteors)):
759 759
760 760 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
761 761 meteorAux = numpy.zeros(16)
762 762
763 763 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
764 764 mHeight = listMeteors[i][0]
765 765 mStart = listMeteors[i][1]
766 766 mPeak = listMeteors[i][2]
767 767 mEnd = listMeteors[i][3]
768 768
769 769 #get the volt data between the start and end times of the meteor
770 770 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
771 771 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
772 772
773 773 #3.6. Phase Difference estimation
774 774 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
775 775
776 776 #3.7. Phase difference removal & meteor start, peak and end times reestimated
777 777 #meteorVolts0.- all Channels, all Profiles
778 778 meteorVolts0 = volts[:,:,mHeight]
779 779 meteorThresh = noise[:,mHeight]*thresholdNoise
780 780 meteorNoise = noise[:,mHeight]
781 781 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
782 782 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
783 783
784 784 #Times reestimation
785 785 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
786 786 if mStart1.size > 0:
787 787 mStart1 = mStart1[-1] + 1
788 788
789 789 else:
790 790 mStart1 = mPeak
791 791
792 792 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
793 793 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
794 794 if mEndDecayTime1.size == 0:
795 795 mEndDecayTime1 = powerNet0.size
796 796 else:
797 797 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
798 798 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
799 799
800 800 #meteorVolts1.- all Channels, from start to end
801 801 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
802 802 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
803 803 if meteorVolts2.shape[1] == 0:
804 804 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
805 805 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
806 806 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
807 807 ##################### END PARAMETERS REESTIMATION #########################
808 808
809 809 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
810 810 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
811 811 if meteorVolts2.shape[1] > 0:
812 812 #Phase Difference re-estimation
813 813 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
814 814 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
815 815 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
816 816 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
817 817 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
818 818
819 819 #Phase Difference RMS
820 820 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
821 821 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
822 822 #Data from Meteor
823 823 mPeak1 = powerNet1.argmax() + mStart1
824 824 mPeakPower1 = powerNet1.max()
825 825 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
826 826 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
827 827 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
828 828 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
829 829 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
830 830 #Vectorize
831 831 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
832 832 meteorAux[7:11] = phaseDiffint[0:4]
833 833
834 834 #Rejection Criterions
835 835 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
836 836 meteorAux[-1] = 17
837 837 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
838 838 meteorAux[-1] = 1
839 839
840 840
841 841 else:
842 842 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
843 843 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
844 844 PowerSeries = 0
845 845
846 846 listMeteors1.append(meteorAux)
847 847 listPowerSeries.append(PowerSeries)
848 848 listVoltageSeries.append(meteorVolts1)
849 849
850 850 return listMeteors1, listPowerSeries, listVoltageSeries
851 851
852 852 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
853 853
854 854 threshError = 10
855 855 #Depending if it is 30 or 50 MHz
856 856 if frequency == 30e6:
857 857 timeLag = 45*10**-3
858 858 else:
859 859 timeLag = 15*10**-3
860 860 lag = numpy.ceil(timeLag/timeInterval)
861 861
862 862 listMeteors1 = []
863 863
864 864 for i in range(len(listMeteors)):
865 865 meteorPower = listPower[i]
866 866 meteorAux = listMeteors[i]
867 867
868 868 if meteorAux[-1] == 0:
869 869
870 870 try:
871 871 indmax = meteorPower.argmax()
872 872 indlag = indmax + lag
873 873
874 874 y = meteorPower[indlag:]
875 875 x = numpy.arange(0, y.size)*timeLag
876 876
877 877 #first guess
878 878 a = y[0]
879 879 tau = timeLag
880 880 #exponential fit
881 881 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
882 882 y1 = self.__exponential_function(x, *popt)
883 883 #error estimation
884 884 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
885 885
886 886 decayTime = popt[1]
887 887 riseTime = indmax*timeInterval
888 888 meteorAux[11:13] = [decayTime, error]
889 889
890 890 #Table items 7, 8 and 11
891 891 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
892 892 meteorAux[-1] = 7
893 893 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
894 894 meteorAux[-1] = 8
895 895 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
896 896 meteorAux[-1] = 11
897 897
898 898
899 899 except:
900 900 meteorAux[-1] = 11
901 901
902 902
903 903 listMeteors1.append(meteorAux)
904 904
905 905 return listMeteors1
906 906
907 907 #Exponential Function
908 908
909 909 def __exponential_function(self, x, a, tau):
910 910 y = a*numpy.exp(-x/tau)
911 911 return y
912 912
913 913 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
914 914
915 915 pairslist1 = list(pairslist)
916 916 pairslist1.append((0,1))
917 917 pairslist1.append((3,4))
918 918 numPairs = len(pairslist1)
919 919 #Time Lag
920 920 timeLag = 45*10**-3
921 921 c = 3e8
922 922 lag = numpy.ceil(timeLag/timeInterval)
923 923 freq = 30e6
924 924
925 925 listMeteors1 = []
926 926
927 927 for i in range(len(listMeteors)):
928 928 meteor = listMeteors[i]
929 929 meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
930 930 if meteor[-1] == 0:
931 931 mStart = listMeteors[i][1]
932 932 mPeak = listMeteors[i][2]
933 933 mLag = mPeak - mStart + lag
934 934
935 935 #get the volt data between the start and end times of the meteor
936 936 meteorVolts = listVolts[i]
937 937 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
938 938
939 939 #Get CCF
940 940 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
941 941
942 942 #Method 2
943 943 slopes = numpy.zeros(numPairs)
944 944 time = numpy.array([-2,-1,1,2])*timeInterval
945 945 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
946 946
947 947 #Correct phases
948 948 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
949 949 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
950 950
951 951 if indDer[0].shape[0] > 0:
952 952 for i in range(indDer[0].shape[0]):
953 953 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
954 954 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
955 955
956 956 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
957 957 for j in range(numPairs):
958 958 fit = stats.linregress(time, angAllCCF[j,:])
959 959 slopes[j] = fit[0]
960 960
961 961 #Remove Outlier
962 962 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
963 963 # slopes = numpy.delete(slopes,indOut)
964 964 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
965 965 # slopes = numpy.delete(slopes,indOut)
966 966
967 967 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
968 968 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
969 969 meteorAux[-2] = radialError
970 970 meteorAux[-3] = radialVelocity
971 971
972 972 #Setting Error
973 973 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
974 974 if numpy.abs(radialVelocity) > 200:
975 975 meteorAux[-1] = 15
976 976 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
977 977 elif radialError > radialStdThresh:
978 978 meteorAux[-1] = 12
979 979
980 980 listMeteors1.append(meteorAux)
981 981 return listMeteors1
982 982
983 983 def __setNewArrays(self, listMeteors, date, heiRang):
984 984
985 985 #New arrays
986 986 arrayMeteors = numpy.array(listMeteors)
987 987 arrayParameters = numpy.zeros((len(listMeteors),10))
988 988
989 989 #Date inclusion
990 990 date = re.findall(r'\((.*?)\)', date)
991 991 date = date[0].split(',')
992 992 date = map(int, date)
993 993 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
994 994 arrayDate = numpy.tile(date, (len(listMeteors), 1))
995 995
996 996 #Meteor array
997 997 arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
998 998 arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
999 999
1000 1000 #Parameters Array
1001 1001 arrayParameters[:,0:3] = arrayMeteors[:,0:3]
1002 1002 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1003 1003
1004 1004 return arrayMeteors, arrayParameters
1005 1005
1006 1006 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1007 1007
1008 1008 arrayAOA = numpy.zeros((phases.shape[0],3))
1009 1009 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1010 1010
1011 1011 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1012 1012 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1013 1013 arrayAOA[:,2] = cosDirError
1014 1014
1015 1015 azimuthAngle = arrayAOA[:,0]
1016 1016 zenithAngle = arrayAOA[:,1]
1017 1017
1018 1018 #Setting Error
1019 1019 #Number 3: AOA not fesible
1020 1020 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1021 1021 error[indInvalid] = 3
1022 1022 #Number 4: Large difference in AOAs obtained from different antenna baselines
1023 1023 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1024 1024 error[indInvalid] = 4
1025 1025 return arrayAOA, error
1026 1026
1027 1027 def __getDirectionCosines(self, arrayPhase, pairsList):
1028 1028
1029 1029 #Initializing some variables
1030 1030 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1031 1031 ang_aux = ang_aux.reshape(1,ang_aux.size)
1032 1032
1033 1033 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1034 1034 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1035 1035
1036 1036
1037 1037 for i in range(2):
1038 1038 #First Estimation
1039 1039 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1040 1040 #Dealias
1041 1041 indcsi = numpy.where(phi0_aux > numpy.pi)
1042 1042 phi0_aux[indcsi] -= 2*numpy.pi
1043 1043 indcsi = numpy.where(phi0_aux < -numpy.pi)
1044 1044 phi0_aux[indcsi] += 2*numpy.pi
1045 1045 #Direction Cosine 0
1046 1046 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1047 1047
1048 1048 #Most-Accurate Second Estimation
1049 1049 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1050 1050 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1051 1051 #Direction Cosine 1
1052 1052 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1053 1053
1054 1054 #Searching the correct Direction Cosine
1055 1055 cosdir0_aux = cosdir0[:,i]
1056 1056 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1057 1057 #Minimum Distance
1058 1058 cosDiff = (cosdir1 - cosdir0_aux)**2
1059 1059 indcos = cosDiff.argmin(axis = 1)
1060 1060 #Saving Value obtained
1061 1061 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1062 1062
1063 1063 return cosdir0, cosdir
1064 1064
1065 1065 def __calculateAOA(self, cosdir, azimuth):
1066 1066 cosdirX = cosdir[:,0]
1067 1067 cosdirY = cosdir[:,1]
1068 1068
1069 1069 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1070 1070 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1071 1071 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1072 1072
1073 1073 return angles
1074 1074
1075 1075 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1076 1076
1077 1077 Ramb = 375 #Ramb = c/(2*PRF)
1078 1078 Re = 6371 #Earth Radius
1079 1079 heights = numpy.zeros(Ranges.shape)
1080 1080
1081 1081 R_aux = numpy.array([0,1,2])*Ramb
1082 1082 R_aux = R_aux.reshape(1,R_aux.size)
1083 1083
1084 1084 Ranges = Ranges.reshape(Ranges.size,1)
1085 1085
1086 1086 Ri = Ranges + R_aux
1087 1087 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1088 1088
1089 1089 #Check if there is a height between 70 and 110 km
1090 1090 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1091 1091 ind_h = numpy.where(h_bool == 1)[0]
1092 1092
1093 1093 hCorr = hi[ind_h, :]
1094 1094 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1095 1095
1096 1096 hCorr = hi[ind_hCorr]
1097 1097 heights[ind_h] = hCorr
1098 1098
1099 1099 #Setting Error
1100 1100 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1101 1101 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1102 1102
1103 1103 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1104 1104 error[indInvalid2] = 14
1105 1105 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1106 1106 error[indInvalid1] = 13
1107 1107
1108 1108 return heights, error
1109 1109
1110 1110
1111 1111 class WindProfiler(Operation):
1112 1112
1113 1113 __isConfig = False
1114 1114
1115 1115 __initime = None
1116 1116 __lastdatatime = None
1117 1117 __integrationtime = None
1118 1118
1119 1119 __buffer = None
1120 1120
1121 1121 __dataReady = False
1122 1122
1123 1123 __firstdata = None
1124 1124
1125 1125 n = None
1126 1126
1127 1127 def __init__(self):
1128 1128 Operation.__init__(self)
1129 1129
1130 def __calculateCosDir(self, elev, azim):
1131 zen = (90 - elev)*numpy.pi/180
1132 azim = azim*numpy.pi/180
1133 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1134 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1135
1136 signX = numpy.sign(numpy.cos(azim))
1137 signY = numpy.sign(numpy.sin(azim))
1138
1139 cosDirX = numpy.copysign(cosDirX, signX)
1140 cosDirY = numpy.copysign(cosDirY, signY)
1141 return cosDirX, cosDirY
1142
1130 1143 def __calculateAngles(self, theta_x, theta_y, azimuth):
1131 1144
1132 1145 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1133 1146 zenith_arr = numpy.arccos(dir_cosw)
1134 1147 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1135 1148
1136 1149 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1137 1150 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1138 1151
1139 1152 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1140 1153
1141 1154 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1142 1155
1143 1156 #
1144 1157 if horOnly:
1145 1158 A = numpy.c_[dir_cosu,dir_cosv]
1146 1159 else:
1147 1160 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1148 1161 A = numpy.asmatrix(A)
1149 1162 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1150 1163
1151 1164 return A1
1152 1165
1153 1166 def __correctValues(self, heiRang, phi, velRadial, SNR):
1154 1167 listPhi = phi.tolist()
1155 1168 maxid = listPhi.index(max(listPhi))
1156 1169 minid = listPhi.index(min(listPhi))
1157 1170
1158 1171 rango = range(len(phi))
1159 1172 # rango = numpy.delete(rango,maxid)
1160 1173
1161 1174 heiRang1 = heiRang*math.cos(phi[maxid])
1162 1175 heiRangAux = heiRang*math.cos(phi[minid])
1163 1176 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1164 1177 heiRang1 = numpy.delete(heiRang1,indOut)
1165 1178
1166 1179 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1167 1180 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1168 1181
1169 1182 for i in rango:
1170 1183 x = heiRang*math.cos(phi[i])
1171 1184 y1 = velRadial[i,:]
1172 1185 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1173 1186
1174 1187 x1 = heiRang1
1175 1188 y11 = f1(x1)
1176 1189
1177 1190 y2 = SNR[i,:]
1178 1191 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1179 1192 y21 = f2(x1)
1180 1193
1181 1194 velRadial1[i,:] = y11
1182 1195 SNR1[i,:] = y21
1183 1196
1184 1197 return heiRang1, velRadial1, SNR1
1185 1198
1186 1199 def __calculateVelUVW(self, A, velRadial):
1187 1200
1188 1201 #Operacion Matricial
1189 1202 # velUVW = numpy.zeros((velRadial.shape[1],3))
1190 1203 # for ind in range(velRadial.shape[1]):
1191 1204 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1192 1205 # velUVW = velUVW.transpose()
1193 1206 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1194 1207 velUVW[:,:] = numpy.dot(A,velRadial)
1195 1208
1196 1209
1197 1210 return velUVW
1198 1211
1199 1212 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1200 1213 """
1201 1214 Function that implements Doppler Beam Swinging (DBS) technique.
1202 1215
1203 1216 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1204 1217 Direction correction (if necessary), Ranges and SNR
1205 1218
1206 1219 Output: Winds estimation (Zonal, Meridional and Vertical)
1207 1220
1208 1221 Parameters affected: Winds, height range, SNR
1209 1222 """
1210 1223 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1211 1224 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1212 1225 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1213 1226
1214 1227 #Calculo de Componentes de la velocidad con DBS
1215 1228 winds = self.__calculateVelUVW(A,velRadial1)
1216 1229
1217 1230 return winds, heiRang1, SNR1
1218 1231
1219 1232 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1220 1233
1221 1234 posx = numpy.asarray(posx)
1222 1235 posy = numpy.asarray(posy)
1223 1236
1224 1237 #Rotacion Inversa para alinear con el azimuth
1225 1238 if azimuth!= None:
1226 1239 azimuth = azimuth*math.pi/180
1227 1240 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1228 1241 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1229 1242 else:
1230 1243 posx1 = posx
1231 1244 posy1 = posy
1232 1245
1233 1246 #Calculo de Distancias
1234 1247 distx = numpy.zeros(pairsCrossCorr.size)
1235 1248 disty = numpy.zeros(pairsCrossCorr.size)
1236 1249 dist = numpy.zeros(pairsCrossCorr.size)
1237 1250 ang = numpy.zeros(pairsCrossCorr.size)
1238 1251
1239 1252 for i in range(pairsCrossCorr.size):
1240 1253 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1241 1254 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1242 1255 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1243 1256 ang[i] = numpy.arctan2(disty[i],distx[i])
1244 1257 #Calculo de Matrices
1245 1258 nPairs = len(pairs)
1246 1259 ang1 = numpy.zeros((nPairs, 2, 1))
1247 1260 dist1 = numpy.zeros((nPairs, 2, 1))
1248 1261
1249 1262 for j in range(nPairs):
1250 1263 dist1[j,0,0] = dist[pairs[j][0]]
1251 1264 dist1[j,1,0] = dist[pairs[j][1]]
1252 1265 ang1[j,0,0] = ang[pairs[j][0]]
1253 1266 ang1[j,1,0] = ang[pairs[j][1]]
1254 1267
1255 1268 return distx,disty, dist1,ang1
1256 1269
1257 1270 def __calculateVelVer(self, phase, lagTRange, _lambda):
1258 1271
1259 1272 Ts = lagTRange[1] - lagTRange[0]
1260 1273 velW = -_lambda*phase/(4*math.pi*Ts)
1261 1274
1262 1275 return velW
1263 1276
1264 1277 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1265 1278 nPairs = tau1.shape[0]
1266 1279 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1267 1280
1268 1281 angCos = numpy.cos(ang)
1269 1282 angSin = numpy.sin(ang)
1270 1283
1271 1284 vel0 = dist*tau1/(2*tau2**2)
1272 1285 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1273 1286 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1274 1287
1275 1288 ind = numpy.where(numpy.isinf(vel))
1276 1289 vel[ind] = numpy.nan
1277 1290
1278 1291 return vel
1279 1292
1280 1293 def __getPairsAutoCorr(self, pairsList, nChannels):
1281 1294
1282 1295 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1283 1296
1284 1297 for l in range(len(pairsList)):
1285 1298 firstChannel = pairsList[l][0]
1286 1299 secondChannel = pairsList[l][1]
1287 1300
1288 1301 #Obteniendo pares de Autocorrelacion
1289 1302 if firstChannel == secondChannel:
1290 1303 pairsAutoCorr[firstChannel] = int(l)
1291 1304
1292 1305 pairsAutoCorr = pairsAutoCorr.astype(int)
1293 1306
1294 1307 pairsCrossCorr = range(len(pairsList))
1295 1308 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1296 1309
1297 1310 return pairsAutoCorr, pairsCrossCorr
1298 1311
1299 1312 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1300 1313 """
1301 1314 Function that implements Spaced Antenna (SA) technique.
1302 1315
1303 1316 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1304 1317 Direction correction (if necessary), Ranges and SNR
1305 1318
1306 1319 Output: Winds estimation (Zonal, Meridional and Vertical)
1307 1320
1308 1321 Parameters affected: Winds
1309 1322 """
1310 1323 #Cross Correlation pairs obtained
1311 1324 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1312 1325 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1313 1326 pairsSelArray = numpy.array(pairsSelected)
1314 1327 pairs = []
1315 1328
1316 1329 #Wind estimation pairs obtained
1317 1330 for i in range(pairsSelArray.shape[0]/2):
1318 1331 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1319 1332 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1320 1333 pairs.append((ind1,ind2))
1321 1334
1322 1335 indtau = tau.shape[0]/2
1323 1336 tau1 = tau[:indtau,:]
1324 1337 tau2 = tau[indtau:-1,:]
1325 1338 tau1 = tau1[pairs,:]
1326 1339 tau2 = tau2[pairs,:]
1327 1340 phase1 = tau[-1,:]
1328 1341
1329 1342 #---------------------------------------------------------------------
1330 1343 #Metodo Directo
1331 1344 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1332 1345 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1333 1346 winds = stats.nanmean(winds, axis=0)
1334 1347 #---------------------------------------------------------------------
1335 1348 #Metodo General
1336 1349 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1337 1350 # #Calculo Coeficientes de Funcion de Correlacion
1338 1351 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1339 1352 # #Calculo de Velocidades
1340 1353 # winds = self.calculateVelUV(F,G,A,B,H)
1341 1354
1342 1355 #---------------------------------------------------------------------
1343 1356 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1344 1357 winds = correctFactor*winds
1345 1358 return winds
1346 1359
1347 1360 def __checkTime(self, currentTime, paramInterval, windsInterval):
1348 1361
1349 1362 dataTime = currentTime + paramInterval
1350 1363 deltaTime = dataTime - self.__initime
1351 1364
1352 1365 if deltaTime >= windsInterval or deltaTime < 0:
1353 1366 self.__dataReady = True
1354 1367 return
1355 1368
1356 1369 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1357 1370 '''
1358 1371 Function that implements winds estimation technique with detected meteors.
1359 1372
1360 1373 Input: Detected meteors, Minimum meteor quantity to wind estimation
1361 1374
1362 1375 Output: Winds estimation (Zonal and Meridional)
1363 1376
1364 1377 Parameters affected: Winds
1365 1378 '''
1366 1379 #Settings
1367 1380 nInt = (heightMax - heightMin)/2
1368 1381 winds = numpy.zeros((2,nInt))*numpy.nan
1369 1382
1370 1383 #Filter errors
1371 1384 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1372 1385 finalMeteor = arrayMeteor[error,:]
1373 1386
1374 1387 #Meteor Histogram
1375 1388 finalHeights = finalMeteor[:,3]
1376 1389 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1377 1390 nMeteorsPerI = hist[0]
1378 1391 heightPerI = hist[1]
1379 1392
1380 1393 #Sort of meteors
1381 1394 indSort = finalHeights.argsort()
1382 1395 finalMeteor2 = finalMeteor[indSort,:]
1383 1396
1384 1397 # Calculating winds
1385 1398 ind1 = 0
1386 1399 ind2 = 0
1387 1400
1388 1401 for i in range(nInt):
1389 1402 nMet = nMeteorsPerI[i]
1390 1403 ind1 = ind2
1391 1404 ind2 = ind1 + nMet
1392 1405
1393 1406 meteorAux = finalMeteor2[ind1:ind2,:]
1394 1407
1395 1408 if meteorAux.shape[0] >= meteorThresh:
1396 1409 vel = meteorAux[:, 7]
1397 1410 zen = meteorAux[:, 5]*numpy.pi/180
1398 1411 azim = meteorAux[:, 4]*numpy.pi/180
1399 1412
1400 1413 n = numpy.cos(zen)
1401 1414 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1402 1415 # l = m*numpy.tan(azim)
1403 1416 l = numpy.sin(zen)*numpy.sin(azim)
1404 1417 m = numpy.sin(zen)*numpy.cos(azim)
1405 1418
1406 1419 A = numpy.vstack((l, m)).transpose()
1407 1420 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1408 1421 windsAux = numpy.dot(A1, vel)
1409 1422
1410 1423 winds[0,i] = windsAux[0]
1411 1424 winds[1,i] = windsAux[1]
1412 1425
1413 1426 return winds, heightPerI[:-1]
1414 1427
1415 1428 def run(self, dataOut, technique, **kwargs):
1416 1429
1417 1430 param = dataOut.data_param
1418 1431 if dataOut.abscissaRange != None:
1419 1432 absc = dataOut.abscissaRange[:-1]
1420 1433 noise = dataOut.noise
1421 1434 heightRange = dataOut.getHeiRange()
1422 1435 SNR = dataOut.SNR
1423 1436
1424 1437 if technique == 'DBS':
1425 1438
1426 theta_x = numpy.array(kwargs['dirCosx'])
1427 theta_y = numpy.array(kwargs['dirCosy'])
1428 azimuth = kwargs['azimuth']
1439 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1440 theta_x = numpy.array(kwargs['dirCosx'])
1441 theta_y = numpy.array(kwargs['dirCosy'])
1442 else:
1443 elev = numpy.array(kwargs['elevation'])
1444 azim = numpy.array(kwargs['azimuth'])
1445 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1446 azimuth = kwargs['correctAzimuth']
1429 1447 if kwargs.has_key('horizontalOnly'):
1430 1448 horizontalOnly = kwargs['horizontalOnly']
1431 1449 else: horizontalOnly = False
1432 1450 if kwargs.has_key('correctFactor'):
1433 1451 correctFactor = kwargs['correctFactor']
1434 1452 else: correctFactor = 1
1435 1453 if kwargs.has_key('channelList'):
1436 1454 channelList = kwargs['channelList']
1437 1455 if len(channelList) == 2:
1438 1456 horizontalOnly = True
1439 1457 arrayChannel = numpy.array(channelList)
1440 1458 param = param[arrayChannel,:,:]
1441 1459 theta_x = theta_x[arrayChannel]
1442 1460 theta_y = theta_y[arrayChannel]
1443 1461
1444 1462 velRadial0 = param[:,1,:] #Radial velocity
1445 1463 dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1446 1464 dataOut.initUtcTime = dataOut.ltctime
1447 1465 dataOut.windsInterval = dataOut.timeInterval
1448 1466
1449 1467 elif technique == 'SA':
1450 1468
1451 1469 #Parameters
1452 1470 position_x = kwargs['positionX']
1453 1471 position_y = kwargs['positionY']
1454 1472 azimuth = kwargs['azimuth']
1455 1473
1456 1474 if kwargs.has_key('crosspairsList'):
1457 1475 pairs = kwargs['crosspairsList']
1458 1476 else:
1459 1477 pairs = None
1460 1478
1461 1479 if kwargs.has_key('correctFactor'):
1462 1480 correctFactor = kwargs['correctFactor']
1463 1481 else:
1464 1482 correctFactor = 1
1465 1483
1466 1484 tau = dataOut.data_param
1467 1485 _lambda = dataOut.C/dataOut.frequency
1468 1486 pairsList = dataOut.pairsList
1469 1487 nChannels = dataOut.nChannels
1470 1488
1471 1489 dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1472 1490 dataOut.initUtcTime = dataOut.ltctime
1473 1491 dataOut.windsInterval = dataOut.timeInterval
1474 1492
1475 1493 elif technique == 'Meteors':
1476 1494 dataOut.flagNoData = True
1477 1495 self.__dataReady = False
1478 1496
1479 1497 if kwargs.has_key('nHours'):
1480 1498 nHours = kwargs['nHours']
1481 1499 else:
1482 1500 nHours = 1
1483 1501
1484 1502 if kwargs.has_key('meteorsPerBin'):
1485 1503 meteorThresh = kwargs['meteorsPerBin']
1486 1504 else:
1487 1505 meteorThresh = 6
1488 1506
1489 1507 if kwargs.has_key('hmin'):
1490 1508 hmin = kwargs['hmin']
1491 1509 else: hmin = 70
1492 1510 if kwargs.has_key('hmax'):
1493 1511 hmax = kwargs['hmax']
1494 1512 else: hmax = 110
1495 1513
1496 1514 dataOut.windsInterval = nHours*3600
1497 1515
1498 1516 if self.__isConfig == False:
1499 1517 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1500 1518 #Get Initial LTC time
1501 1519 self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1502 1520 self.__isConfig = True
1503 1521
1504 1522 if self.__buffer == None:
1505 1523 self.__buffer = dataOut.data_param
1506 1524 self.__firstdata = copy.copy(dataOut)
1507 1525
1508 1526 else:
1509 1527 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1510 1528
1511 1529 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
1512 1530
1513 1531 if self.__dataReady:
1514 1532 dataOut.initUtcTime = self.__initime
1515 1533 self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
1516 1534
1517 1535 dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1518 1536 dataOut.flagNoData = False
1519 1537 self.__buffer = None
1520 1538
1521 1539 return No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,144 +1,144
1 1 # DIAS 19 Y 20 FEB 2014
2 2 # Comprobacion de Resultados DBS con SA
3 3
4 4 import os, sys
5 5
6 6 path = os.path.split(os.getcwd())[0]
7 7 sys.path.append(path)
8 8
9 9 from controller import *
10 10
11 11 desc = "DBS Experiment Test"
12 12 filename = "DBStest.xml"
13 13
14 14 controllerObj = Project()
15 15
16 16 controllerObj.setup(id = '191', name='test01', description=desc)
17 17
18 18 #Experimentos
19 19
20 20 #2014050 19 Feb 2014
21 21 # path = '/home/soporte/Documents/MST_Data/DBS/d2014050'
22 22 # pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/'
23 23 # xmin = '15.5'
24 24 # xmax = '23.99999999'
25 25 # startTime = '17:25:00'
26 26 # filehdf5 = "DBS_2014050.hdf5"
27 27
28 28 #2014051 20 Feb 2014
29 29 path = '/home/soporte/Data/MST/DBS/d2014051'
30 30 pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/'
31 31 xmin = '0.0'
32 32 xmax = '8.0'
33 33 startTime = '00:00:00'
34 34 filehdf5 = "DBS_2014051.hdf5"
35 35
36 36
37 37
38 38 #------------------------------------------------------------------------------------------------
39 39 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
40 40 path=path,
41 41 startDate='2014/01/31',
42 42 endDate='2014/03/31',
43 43 startTime=startTime,
44 44 endTime='23:59:59',
45 45 online=0,
46 46 delay=5,
47 47 walk=0)
48 48
49 49 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
50 50
51 51
52 52 #--------------------------------------------------------------------------------------------------
53 53
54 54 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
55 55
56 56 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
57 57
58 58 opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
59 59 opObj11.addParameter(name='n', value='256', format='int')
60 60 # opObj11.addParameter(name='n', value='16', format='int')
61 61
62 62 opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
63 63 opObj11.addParameter(name='minIndex', value='10', format='float')
64 64 opObj11.addParameter(name='maxIndex', value='60', format='float')
65 65
66 66 #---------------------------------------------------------------------------------------------------
67 67
68 68 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
69 69 procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int')
70 70 procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int')
71 71 # procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int')
72 72 procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
73 73
74 74 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
75 75 opObj11.addParameter(name='n', value='5', format='int')
76 76
77 77 opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
78 78 opObj14.addParameter(name='id', value='1', format='int')
79 79 opObj14.addParameter(name='wintitle', value='Con interf', format='str')
80 80 opObj14.addParameter(name='save', value='1', format='bool')
81 81 opObj14.addParameter(name='figpath', value=pathFigure, format='str')
82 82 opObj14.addParameter(name='zmin', value='5', format='int')
83 83 opObj14.addParameter(name='zmax', value='90', format='int')
84 84
85 85 opObj12 = procUnitConfObj1.addOperation(name='removeInterference')
86 86 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
87 87 opObj13.addParameter(name='mode', value='1', format='int')
88 88
89 89 opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
90 90 opObj12.addParameter(name='id', value='2', format='int')
91 91 opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
92 92 opObj12.addParameter(name='save', value='1', format='bool')
93 93 opObj12.addParameter(name='figpath', value = pathFigure, format='str')
94 94 opObj12.addParameter(name='xmin', value=xmin, format='float')
95 95 opObj12.addParameter(name='xmax', value=xmax, format='float')
96 96 opObj12.addParameter(name='zmin', value='5', format='int')
97 97 opObj12.addParameter(name='zmax', value='90', format='int')
98 98
99 99 #--------------------------------------------------------------------------------------------------
100 100
101 101 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
102 102 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
103 103
104 104 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
105 105 opObj21.addParameter(name='id', value='3', format='int')
106 106 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
107 107 opObj21.addParameter(name='save', value='1', format='bool')
108 108 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
109 109 opObj21.addParameter(name='zmin', value='5', format='int')
110 110 opObj21.addParameter(name='zmax', value='90', format='int')
111 111
112 112 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
113 113 opObj22.addParameter(name='technique', value='DBS', format='str')
114 opObj22.addParameter(name='azimuth', value='51.06', format='float')
114 opObj22.addParameter(name='correctAzimuth', value='51.06', format='float')
115 115 opObj22.addParameter(name='correctFactor', value='-1', format='float')
116 116 opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist')
117 117 opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist')
118 118 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
119 119 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
120 120
121 121 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
122 122 opObj23.addParameter(name='id', value='4', format='int')
123 123 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
124 124 opObj23.addParameter(name='save', value='1', format='bool')
125 125 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
126 126 opObj23.addParameter(name='zmin', value='-10', format='int')
127 127 opObj23.addParameter(name='zmax', value='10', format='int')
128 128 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
129 129 opObj23.addParameter(name='zmax_ver', value='80', format='float')
130 130 opObj23.addParameter(name='SNRmin', value='-10', format='int')
131 131 opObj23.addParameter(name='SNRmax', value='60', format='int')
132 132 opObj23.addParameter(name='SNRthresh', value='0', format='float')
133 133 opObj23.addParameter(name='xmin', value=xmin, format='float')
134 134 opObj23.addParameter(name='xmax', value=xmax, format='float')
135 135
136 136 #--------------------------------------------------------------------------------------------------
137 137 print "Escribiendo el archivo XML"
138 138 controllerObj.writeXml(filename)
139 139 print "Leyendo el archivo XML"
140 140 controllerObj.readXml(filename)
141 141
142 142 controllerObj.createObjects()
143 143 controllerObj.connectObjects()
144 144 controllerObj.run() No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,131 +1,146
1 1 import os, sys
2 2
3 3 path = os.path.split(os.getcwd())[0]
4 4 sys.path.append(path)
5 5
6 6 from controller import *
7 7
8 8 desc = "AMISR Experiment"
9 9
10 10 filename = "amisr_reader.xml"
11 11
12 12 controllerObj = Project()
13 13
14 14 controllerObj.setup(id = '191', name='test01', description=desc)
15 15
16 16
17 17 path = os.path.join(os.environ['HOME'],'Development/amisr/data')
18 18 path = '/home/soporte/Data/AMISR'
19 19 figpath = os.path.join(os.environ['HOME'],'Pictures/amisr')
20 20
21 21 pathFigure = '/home/soporte/workspace/Graficos/DBS/amisr/'
22 xmin = '12.0'
23 xmax = '14.0'
22 xmin = '16.0'
23 xmax = '17.0'
24 24
25 25 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader',
26 26 path=path,
27 27 startDate='2014/10/21',
28 28 endDate='2014/10/21',
29 29 startTime='00:00:00',
30 30 endTime='23:59:59',
31 31 walk=1,
32 32 timezone='lt')
33 33
34 34 #AMISR Processing Unit
35 35 procUnitAMISRBeam0 = controllerObj.addProcUnit(datatype='AMISRProc', inputId=readUnitConfObj.getId())
36 36
37 37 opObj11 = procUnitAMISRBeam0.addOperation(name='PrintInfo', optype='other')
38 38
39 39 #Reshaper
40 40 opObj11 = procUnitAMISRBeam0.addOperation(name='ProfileToChannels', optype='other')
41 41
42 42
43 43 #Beam Selector
44 44 #opObj11 = procUnitAMISRBeam0.addOperation(name='BeamSelector', optype='other')
45 45 #opObj11.addParameter(name='beam', value='0', format='int')
46 46
47 47 #Voltage Processing Unit
48 48 procUnitConfObjBeam0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=procUnitAMISRBeam0.getId())
49 opObj10 = procUnitConfObjBeam0.addOperation(name='setRadarFrequency')
50 opObj10.addParameter(name='frequency', value='445e6', format='float')
51
52 opObj12 = procUnitConfObjBeam0.addOperation(name='selectHeights')
53 opObj12.addParameter(name='minHei', value='0', format='float')
54 opObj12.addParameter(name='maxHei', value='10', format='float')
49 55 #Coherent Integration
50 56 opObj11 = procUnitConfObjBeam0.addOperation(name='CohInt', optype='other')
51 57 opObj11.addParameter(name='n', value='8', format='int')
52 58 #Spectra Unit Processing, getting spectras with nProfiles and nFFTPoints
53 59 procUnitConfObjSpectraBeam0 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjBeam0.getId())
54 60 procUnitConfObjSpectraBeam0.addParameter(name='nFFTPoints', value=32, format='int')
55 61 procUnitConfObjSpectraBeam0.addParameter(name='nProfiles', value=32, format='int')
62
63 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='IncohInt', optype='other')
64 opObj11.addParameter(name='n', value='16', format='int')
65
56 66 #RemoveDc
57 67 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='removeDC')
58 68
59 69 #Noise Estimation
60 70 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='getNoise')
61 71 opObj11.addParameter(name='minHei', value='5', format='float')
62 72 opObj11.addParameter(name='maxHei', value='20', format='float')
63 73
64 74 #SpectraPlot
65 75 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='SpectraPlot', optype='other')
66 76 opObj11.addParameter(name='id', value='100', format='int')
67 77 opObj11.addParameter(name='wintitle', value='AMISR Beam 0', format='str')
68 78 opObj11.addParameter(name='zmin', value='30', format='int')
69 79 opObj11.addParameter(name='zmax', value='80', format='int')
70
80 opObj11.addParameter(name='save', value='1', format='bool')
81 opObj11.addParameter(name='figpath', value = pathFigure, format='str')
71 82 #RTIPlot
72 83 #title0 = 'RTI AMISR Beam 0'
73 #opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other')
74 #opObj11.addParameter(name='id', value='200', format='int')
75 #opObj11.addParameter(name='wintitle', value=title0, format='str')
76 #opObj11.addParameter(name='showprofile', value='0', format='int')
77 ##Setting RTI time using xmin,xmax
78 #opObj11.addParameter(name='xmin', value='15', format='int')
79 #opObj11.addParameter(name='xmax', value='23', format='int')
80 #Setting dB range with zmin, zmax
81 #opObj11.addParameter(name='zmin', value='45', format='int')
82 #opObj11.addParameter(name='zmax', value='70', format='int')
83 #Save RTI
84 #figfile0 = 'amisr_rti_beam0.png'
85 #opObj11.addParameter(name='figpath', value=figpath, format='str')
86 #opObj11.addParameter(name='figfile', value=figfile0, format='str')
84 opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other')
85 opObj11.addParameter(name='id', value='200', format='int')
86 # opObj11.addParameter(name='wintitle', value=title0, format='str')
87 opObj11.addParameter(name='showprofile', value='0', format='int')
88 #Setting RTI time using xmin,xmax
89 opObj11.addParameter(name='xmin', value=xmin, format='float')
90 opObj11.addParameter(name='xmax', value=xmax, format='float')
91 # Setting dB range with zmin, zmax
92 opObj11.addParameter(name='zmin', value='45', format='int')
93 opObj11.addParameter(name='zmax', value='70', format='int')
94 # Save RTI
95 # figfile0 = 'amisr_rti_beam0.png'
96 # opObj11.addParameter(name='figpath', value=figpath, format='str')
97 # opObj11.addParameter(name='figfile', value=figfile0, format='str')
87 98
88 99 #-----------------------------------------------------------------------------------------------
89 100 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjSpectraBeam0 .getId())
90 101 opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
91 102
92 # opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
93 # opObj21.addParameter(name='id', value='3', format='int')
94 # opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
95 # opObj21.addParameter(name='save', value='1', format='bool')
96 # opObj21.addParameter(name='figpath', value=pathFigure, format='str')
97 # opObj21.addParameter(name='zmin', value='5', format='int')
98 # opObj21.addParameter(name='zmax', value='90', format='int')
103 opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
104 opObj21.addParameter(name='id', value='3', format='int')
105 opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
106 opObj21.addParameter(name='save', value='1', format='bool')
107 opObj21.addParameter(name='figpath', value=pathFigure, format='str')
108 opObj21.addParameter(name='zmin', value='30', format='int')
109 opObj21.addParameter(name='zmax', value='80', format='int')
99 110
100 111 opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
101 112 opObj22.addParameter(name='technique', value='DBS', format='str')
102 opObj22.addParameter(name='azimuth', value='51.06', format='float')
113 opObj22.addParameter(name='correctAzimuth', value='51.06', format='float')
103 114 opObj22.addParameter(name='correctFactor', value='-1', format='float')
104 opObj22.addParameter(name='dirCosx', value='9.63770247e-01, 5.93066547e-17, 1, 5.93066547e-17,-9.68583161e-01', format='floatlist')
105 opObj22.addParameter(name='dirCosy', value=' 0,-0.96858316,0,0.96858316,0', format='floatlist')
115 opObj22.addParameter(name='azimuth', value='0,180,90,0,-90', format='floatlist')
116 opObj22.addParameter(name='elevation', value='74.53,75.60,75.60,90.0,75.60', format='floatlist')
106 117 # opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
107 118 # opObj22.addParameter(name='channelList', value='1,2', format='intlist')
108 119
109 120 opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
110 121 opObj23.addParameter(name='id', value='4', format='int')
111 122 opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
112 123 opObj23.addParameter(name='save', value='1', format='bool')
113 124 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
114 opObj23.addParameter(name='zmin', value='-10', format='int')
115 opObj23.addParameter(name='zmax', value='10', format='int')
116 opObj23.addParameter(name='zmin_ver', value='-80', format='float')
117 opObj23.addParameter(name='zmax_ver', value='80', format='float')
125 opObj23.addParameter(name='zmin', value='-20', format='int')
126 opObj23.addParameter(name='zmax', value='20', format='int')
127 opObj23.addParameter(name='zmin_ver', value='-100', format='float')
128 opObj23.addParameter(name='zmax_ver', value='100', format='float')
118 129 opObj23.addParameter(name='SNRmin', value='-10', format='int')
119 130 opObj23.addParameter(name='SNRmax', value='60', format='int')
120 131 opObj23.addParameter(name='SNRthresh', value='0', format='float')
121 132 opObj23.addParameter(name='xmin', value=xmin, format='float')
122 133 opObj23.addParameter(name='xmax', value=xmax, format='float')
123 134
124 135 print "Escribiendo el archivo XML"
125 136 controllerObj.writeXml(filename)
126 137 print "Leyendo el archivo XML"
127 138 controllerObj.readXml(filename)
128 139
129 140 controllerObj.createObjects()
130 141 controllerObj.connectObjects()
131 142 controllerObj.run()
143
144 #21 3 pm
145
146
General Comments 0
You need to be logged in to leave comments. Login now