##// END OF EJS Templates
Final Commit! GG WP, no bugs.
ebocanegra -
r1158:a0a227f9ae68 claire_proc
parent child
Show More
@@ -1,795 +1,803
1 1 import os, sys
2 2 import glob
3 3 import fnmatch
4 4 import datetime
5 5 import time
6 6 import re
7 7 import h5py
8 8 import numpy
9 9 import matplotlib.pyplot as plt
10 10
11 11 import pylab as plb
12 12 from scipy.optimize import curve_fit
13 13 from scipy import asarray as ar, exp
14 14 from scipy import stats
15 15
16 16 from numpy.ma.core import getdata
17 17
18 18 SPEED_OF_LIGHT = 299792458
19 19 SPEED_OF_LIGHT = 3e8
20 20
21 21 try:
22 22 from gevent import sleep
23 23 except:
24 24 from time import sleep
25 25
26 26 from schainpy.model.data.jrodata import Spectra
27 27 #from schainpy.model.data.BLTRheaderIO import FileHeader, RecordHeader
28 28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
29 29 #from schainpy.model.io.jroIO_bltr import BLTRReader
30 30 from numpy import imag, shape, NaN
31 31
32 32 from jroIO_base import JRODataReader
33 33
34 34
35 35 class Header(object):
36 36
37 37 def __init__(self):
38 38 raise NotImplementedError
39 39
40 40
41 41 def read(self):
42 42
43 43 raise NotImplementedError
44 44
45 45 def write(self):
46 46
47 47 raise NotImplementedError
48 48
49 49 def printInfo(self):
50 50
51 51 message = "#"*50 + "\n"
52 52 message += self.__class__.__name__.upper() + "\n"
53 53 message += "#"*50 + "\n"
54 54
55 55 keyList = self.__dict__.keys()
56 56 keyList.sort()
57 57
58 58 for key in keyList:
59 59 message += "%s = %s" %(key, self.__dict__[key]) + "\n"
60 60
61 61 if "size" not in keyList:
62 62 attr = getattr(self, "size")
63 63
64 64 if attr:
65 65 message += "%s = %s" %("size", attr) + "\n"
66 66
67 67 #print message
68 68
69 69
70 70
71 71
72 72
73 73 FILE_STRUCTURE = numpy.dtype([ #HEADER 48bytes
74 74 ('FileMgcNumber','<u4'), #0x23020100
75 75 ('nFDTdataRecors','<u4'), #No Of FDT data records in this file (0 or more)
76 76 ('OffsetStartHeader','<u4'),
77 77 ('RadarUnitId','<u4'),
78 78 ('SiteName',numpy.str_,32), #Null terminated
79 79 ])
80 80
81 81 class FileHeaderBLTR(Header):
82 82
83 83 def __init__(self):
84 84
85 85 self.FileMgcNumber= 0 #0x23020100
86 86 self.nFDTdataRecors=0 #No Of FDT data records in this file (0 or more)
87 87 self.RadarUnitId= 0
88 88 self.OffsetStartHeader=0
89 89 self.SiteName= ""
90 90 self.size = 48
91 91
92 92 def FHread(self, fp):
93 93 #try:
94 94 startFp = open(fp,"rb")
95 95
96 96 header = numpy.fromfile(startFp, FILE_STRUCTURE,1)
97 97
98 98 print ' '
99 99 print 'puntero file header', startFp.tell()
100 100 print ' '
101 101
102 102
103 103 ''' numpy.fromfile(file, dtype, count, sep='')
104 104 file : file or str
105 105 Open file object or filename.
106 106
107 107 dtype : data-type
108 108 Data type of the returned array. For binary files, it is used to determine
109 109 the size and byte-order of the items in the file.
110 110
111 111 count : int
112 112 Number of items to read. -1 means all items (i.e., the complete file).
113 113
114 114 sep : str
115 115 Separator between items if file is a text file. Empty ("") separator means
116 116 the file should be treated as binary. Spaces (" ") in the separator match zero
117 117 or more whitespace characters. A separator consisting only of spaces must match
118 118 at least one whitespace.
119 119
120 120 '''
121 121
122 122
123 123
124 124 self.FileMgcNumber= hex(header['FileMgcNumber'][0])
125 125 self.nFDTdataRecors=int(header['nFDTdataRecors'][0]) #No Of FDT data records in this file (0 or more)
126 126 self.RadarUnitId= int(header['RadarUnitId'][0])
127 127 self.OffsetStartHeader= int(header['OffsetStartHeader'][0])
128 128 self.SiteName= str(header['SiteName'][0])
129 129
130 130 #print 'Numero de bloques', self.nFDTdataRecors
131 131
132 132
133 133 if self.size <48:
134 134 return 0
135 135
136 136 return 1
137 137
138 138
139 139 def write(self, fp):
140 140
141 141 headerTuple = (self.FileMgcNumber,
142 142 self.nFDTdataRecors,
143 143 self.RadarUnitId,
144 144 self.SiteName,
145 145 self.size)
146 146
147 147
148 148 header = numpy.array(headerTuple, FILE_STRUCTURE)
149 149 # numpy.array(object, dtype=None, copy=True, order=None, subok=False, ndmin=0)
150 150 header.tofile(fp)
151 151 ''' ndarray.tofile(fid, sep, format) Write array to a file as text or binary (default).
152 152
153 153 fid : file or str
154 154 An open file object, or a string containing a filename.
155 155
156 156 sep : str
157 157 Separator between array items for text output. If "" (empty), a binary file is written,
158 158 equivalent to file.write(a.tobytes()).
159 159
160 160 format : str
161 161 Format string for text file output. Each entry in the array is formatted to text by
162 162 first converting it to the closest Python type, and then using "format" % item.
163 163
164 164 '''
165 165
166 166 return 1
167 167
168 168
169 169
170 170
171 171
172 172 RECORD_STRUCTURE = numpy.dtype([ #RECORD HEADER 180+20N bytes
173 173 ('RecMgcNumber','<u4'), #0x23030001
174 174 ('RecCounter','<u4'), #Record counter(0,1, ...)
175 175 ('Off2StartNxtRec','<u4'), #Offset to start of next record form start of this record
176 176 ('Off2StartData','<u4'), #Offset to start of data from start of this record
177 177 ('nUtime','<i4'), #Epoch time stamp of start of acquisition (seconds)
178 178 ('nMilisec','<u4'), #Millisecond component of time stamp (0,...,999)
179 179 ('ExpTagName',numpy.str_,32), #Experiment tag name (null terminated)
180 180 ('ExpComment',numpy.str_,32), #Experiment comment (null terminated)
181 181 ('SiteLatDegrees','<f4'), #Site latitude (from GPS) in degrees (positive implies North)
182 182 ('SiteLongDegrees','<f4'), #Site longitude (from GPS) in degrees (positive implies East)
183 183 ('RTCgpsStatus','<u4'), #RTC GPS engine status (0=SEEK, 1=LOCK, 2=NOT FITTED, 3=UNAVAILABLE)
184 184 ('TransmitFrec','<u4'), #Transmit frequency (Hz)
185 185 ('ReceiveFrec','<u4'), #Receive frequency
186 186 ('FirstOsciFrec','<u4'), #First local oscillator frequency (Hz)
187 187 ('Polarisation','<u4'), #(0="O", 1="E", 2="linear 1", 3="linear2")
188 188 ('ReceiverFiltSett','<u4'), #Receiver filter settings (0,1,2,3)
189 189 ('nModesInUse','<u4'), #Number of modes in use (1 or 2)
190 190 ('DualModeIndex','<u4'), #Dual Mode index number for these data (0 or 1)
191 191 ('DualModeRange','<u4'), #Dual Mode range correction for these data (m)
192 192 ('nDigChannels','<u4'), #Number of digital channels acquired (2*N)
193 193 ('SampResolution','<u4'), #Sampling resolution (meters)
194 194 ('nHeights','<u4'), #Number of range gates sampled
195 195 ('StartRangeSamp','<u4'), #Start range of sampling (meters)
196 196 ('PRFhz','<u4'), #PRF (Hz)
197 197 ('nCohInt','<u4'), #Integrations
198 198 ('nProfiles','<u4'), #Number of data points transformed
199 199 ('nChannels','<u4'), #Number of receive beams stored in file (1 or N)
200 200 ('nIncohInt','<u4'), #Number of spectral averages
201 201 ('FFTwindowingInd','<u4'), #FFT windowing index (0 = no window)
202 202 ('BeamAngleAzim','<f4'), #Beam steer angle (azimuth) in degrees (clockwise from true North)
203 203 ('BeamAngleZen','<f4'), #Beam steer angle (zenith) in degrees (0=> vertical)
204 204 ('AntennaCoord0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
205 205 ('AntennaAngl0','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
206 206 ('AntennaCoord1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
207 207 ('AntennaAngl1','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
208 208 ('AntennaCoord2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
209 209 ('AntennaAngl2','<f4'), #Antenna coordinates (Range(meters), Bearing(degrees)) - N pairs
210 210 ('RecPhaseCalibr0','<f4'), #Receiver phase calibration (degrees) - N values
211 211 ('RecPhaseCalibr1','<f4'), #Receiver phase calibration (degrees) - N values
212 212 ('RecPhaseCalibr2','<f4'), #Receiver phase calibration (degrees) - N values
213 213 ('RecAmpCalibr0','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
214 214 ('RecAmpCalibr1','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
215 215 ('RecAmpCalibr2','<f4'), #Receiver amplitude calibration (ratio relative to receiver one) - N values
216 216 ('ReceiverGaindB0','<i4'), #Receiver gains in dB - N values
217 217 ('ReceiverGaindB1','<i4'), #Receiver gains in dB - N values
218 218 ('ReceiverGaindB2','<i4'), #Receiver gains in dB - N values
219 219 ])
220 220
221 221
222 222 class RecordHeaderBLTR(Header):
223 223
224 224 def __init__(self, RecMgcNumber=None, RecCounter= 0, Off2StartNxtRec= 811248,
225 225 nUtime= 0, nMilisec= 0, ExpTagName= None,
226 226 ExpComment=None, SiteLatDegrees=0, SiteLongDegrees= 0,
227 227 RTCgpsStatus= 0, TransmitFrec= 0, ReceiveFrec= 0,
228 228 FirstOsciFrec= 0, Polarisation= 0, ReceiverFiltSett= 0,
229 229 nModesInUse= 0, DualModeIndex= 0, DualModeRange= 0,
230 230 nDigChannels= 0, SampResolution= 0, nHeights= 0,
231 231 StartRangeSamp= 0, PRFhz= 0, nCohInt= 0,
232 232 nProfiles= 0, nChannels= 0, nIncohInt= 0,
233 233 FFTwindowingInd= 0, BeamAngleAzim= 0, BeamAngleZen= 0,
234 234 AntennaCoord0= 0, AntennaCoord1= 0, AntennaCoord2= 0,
235 235 RecPhaseCalibr0= 0, RecPhaseCalibr1= 0, RecPhaseCalibr2= 0,
236 236 RecAmpCalibr0= 0, RecAmpCalibr1= 0, RecAmpCalibr2= 0,
237 237 AntennaAngl0=0, AntennaAngl1=0, AntennaAngl2=0,
238 238 ReceiverGaindB0= 0, ReceiverGaindB1= 0, ReceiverGaindB2= 0, Off2StartData=0, OffsetStartHeader=0):
239 239
240 240 self.RecMgcNumber = RecMgcNumber #0x23030001
241 241 self.RecCounter = RecCounter
242 242 self.Off2StartNxtRec = Off2StartNxtRec
243 243 self.Off2StartData = Off2StartData
244 244 self.nUtime = nUtime
245 245 self.nMilisec = nMilisec
246 246 self.ExpTagName = ExpTagName
247 247 self.ExpComment = ExpComment
248 248 self.SiteLatDegrees = SiteLatDegrees
249 249 self.SiteLongDegrees = SiteLongDegrees
250 250 self.RTCgpsStatus = RTCgpsStatus
251 251 self.TransmitFrec = TransmitFrec
252 252 self.ReceiveFrec = ReceiveFrec
253 253 self.FirstOsciFrec = FirstOsciFrec
254 254 self.Polarisation = Polarisation
255 255 self.ReceiverFiltSett = ReceiverFiltSett
256 256 self.nModesInUse = nModesInUse
257 257 self.DualModeIndex = DualModeIndex
258 258 self.DualModeRange = DualModeRange
259 259 self.nDigChannels = nDigChannels
260 260 self.SampResolution = SampResolution
261 261 self.nHeights = nHeights
262 262 self.StartRangeSamp = StartRangeSamp
263 263 self.PRFhz = PRFhz
264 264 self.nCohInt = nCohInt
265 265 self.nProfiles = nProfiles
266 266 self.nChannels = nChannels
267 267 self.nIncohInt = nIncohInt
268 268 self.FFTwindowingInd = FFTwindowingInd
269 269 self.BeamAngleAzim = BeamAngleAzim
270 270 self.BeamAngleZen = BeamAngleZen
271 271 self.AntennaCoord0 = AntennaCoord0
272 272 self.AntennaAngl0 = AntennaAngl0
273 273 self.AntennaAngl1 = AntennaAngl1
274 274 self.AntennaAngl2 = AntennaAngl2
275 275 self.AntennaCoord1 = AntennaCoord1
276 276 self.AntennaCoord2 = AntennaCoord2
277 277 self.RecPhaseCalibr0 = RecPhaseCalibr0
278 278 self.RecPhaseCalibr1 = RecPhaseCalibr1
279 279 self.RecPhaseCalibr2 = RecPhaseCalibr2
280 280 self.RecAmpCalibr0 = RecAmpCalibr0
281 281 self.RecAmpCalibr1 = RecAmpCalibr1
282 282 self.RecAmpCalibr2 = RecAmpCalibr2
283 283 self.ReceiverGaindB0 = ReceiverGaindB0
284 284 self.ReceiverGaindB1 = ReceiverGaindB1
285 285 self.ReceiverGaindB2 = ReceiverGaindB2
286 286 self.OffsetStartHeader = 48
287 287
288 288
289 289
290 290 def RHread(self, fp):
291 291 #print fp
292 292 #startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
293 293 startFp = open(fp,"rb") #The method tell() returns the current position of the file read/write pointer within the file.
294 294 #RecCounter=0
295 295 #Off2StartNxtRec=811248
296 296 OffRHeader= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
297 297 print ' '
298 298 print 'puntero Record Header', startFp.tell()
299 299 print ' '
300 300
301 301
302 302 startFp.seek(OffRHeader, os.SEEK_SET)
303 303
304 304 print ' '
305 305 print 'puntero Record Header con seek', startFp.tell()
306 306 print ' '
307 307
308 308 #print 'Posicion del bloque: ',OffRHeader
309 309
310 310 header = numpy.fromfile(startFp,RECORD_STRUCTURE,1)
311 311
312 312 print ' '
313 313 print 'puntero Record Header con seek', startFp.tell()
314 314 print ' '
315 315
316 316 print ' '
317 317 #
318 318 #print 'puntero Record Header despues de seek', header.tell()
319 319 print ' '
320 320
321 321 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) #0x23030001
322 322 self.RecCounter = int(header['RecCounter'][0])
323 323 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
324 324 self.Off2StartData = int(header['Off2StartData'][0])
325 325 self.nUtime = header['nUtime'][0]
326 326 self.nMilisec = header['nMilisec'][0]
327 327 self.ExpTagName = str(header['ExpTagName'][0])
328 328 self.ExpComment = str(header['ExpComment'][0])
329 329 self.SiteLatDegrees = header['SiteLatDegrees'][0]
330 330 self.SiteLongDegrees = header['SiteLongDegrees'][0]
331 331 self.RTCgpsStatus = header['RTCgpsStatus'][0]
332 332 self.TransmitFrec = header['TransmitFrec'][0]
333 333 self.ReceiveFrec = header['ReceiveFrec'][0]
334 334 self.FirstOsciFrec = header['FirstOsciFrec'][0]
335 335 self.Polarisation = header['Polarisation'][0]
336 336 self.ReceiverFiltSett = header['ReceiverFiltSett'][0]
337 337 self.nModesInUse = header['nModesInUse'][0]
338 338 self.DualModeIndex = header['DualModeIndex'][0]
339 339 self.DualModeRange = header['DualModeRange'][0]
340 340 self.nDigChannels = header['nDigChannels'][0]
341 341 self.SampResolution = header['SampResolution'][0]
342 342 self.nHeights = header['nHeights'][0]
343 343 self.StartRangeSamp = header['StartRangeSamp'][0]
344 344 self.PRFhz = header['PRFhz'][0]
345 345 self.nCohInt = header['nCohInt'][0]
346 346 self.nProfiles = header['nProfiles'][0]
347 347 self.nChannels = header['nChannels'][0]
348 348 self.nIncohInt = header['nIncohInt'][0]
349 349 self.FFTwindowingInd = header['FFTwindowingInd'][0]
350 350 self.BeamAngleAzim = header['BeamAngleAzim'][0]
351 351 self.BeamAngleZen = header['BeamAngleZen'][0]
352 352 self.AntennaCoord0 = header['AntennaCoord0'][0]
353 353 self.AntennaAngl0 = header['AntennaAngl0'][0]
354 354 self.AntennaCoord1 = header['AntennaCoord1'][0]
355 355 self.AntennaAngl1 = header['AntennaAngl1'][0]
356 356 self.AntennaCoord2 = header['AntennaCoord2'][0]
357 357 self.AntennaAngl2 = header['AntennaAngl2'][0]
358 358 self.RecPhaseCalibr0 = header['RecPhaseCalibr0'][0]
359 359 self.RecPhaseCalibr1 = header['RecPhaseCalibr1'][0]
360 360 self.RecPhaseCalibr2 = header['RecPhaseCalibr2'][0]
361 361 self.RecAmpCalibr0 = header['RecAmpCalibr0'][0]
362 362 self.RecAmpCalibr1 = header['RecAmpCalibr1'][0]
363 363 self.RecAmpCalibr2 = header['RecAmpCalibr2'][0]
364 364 self.ReceiverGaindB0 = header['ReceiverGaindB0'][0]
365 365 self.ReceiverGaindB1 = header['ReceiverGaindB1'][0]
366 366 self.ReceiverGaindB2 = header['ReceiverGaindB2'][0]
367 367
368 368 self.ipp= 0.5*(SPEED_OF_LIGHT/self.PRFhz)
369 369
370 370 self.RHsize = 180+20*self.nChannels
371 371 self.Datasize= self.nProfiles*self.nChannels*self.nHeights*2*4
372 372 #print 'Datasize',self.Datasize
373 373 endFp = self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec
374 374
375 375 print '=============================================='
376 376 print 'RecMgcNumber ',self.RecMgcNumber
377 377 print 'RecCounter ',self.RecCounter
378 378 print 'Off2StartNxtRec ',self.Off2StartNxtRec
379 379 print 'Off2StartData ',self.Off2StartData
380 380 print 'Range Resolution ',self.SampResolution
381 381 print 'First Height ',self.StartRangeSamp
382 382 print 'PRF (Hz) ',self.PRFhz
383 383 print 'Heights (K) ',self.nHeights
384 384 print 'Channels (N) ',self.nChannels
385 385 print 'Profiles (J) ',self.nProfiles
386 386 print 'iCoh ',self.nCohInt
387 387 print 'iInCoh ',self.nIncohInt
388 388 print 'BeamAngleAzim ',self.BeamAngleAzim
389 389 print 'BeamAngleZen ',self.BeamAngleZen
390 390
391 391 #print 'ModoEnUso ',self.DualModeIndex
392 392 #print 'UtcTime ',self.nUtime
393 393 #print 'MiliSec ',self.nMilisec
394 394 #print 'Exp TagName ',self.ExpTagName
395 395 #print 'Exp Comment ',self.ExpComment
396 396 #print 'FFT Window Index ',self.FFTwindowingInd
397 397 #print 'N Dig. Channels ',self.nDigChannels
398 398 print 'Size de bloque ',self.RHsize
399 399 print 'DataSize ',self.Datasize
400 400 print 'BeamAngleAzim ',self.BeamAngleAzim
401 401 #print 'AntennaCoord0 ',self.AntennaCoord0
402 402 #print 'AntennaAngl0 ',self.AntennaAngl0
403 403 #print 'AntennaCoord1 ',self.AntennaCoord1
404 404 #print 'AntennaAngl1 ',self.AntennaAngl1
405 405 #print 'AntennaCoord2 ',self.AntennaCoord2
406 406 #print 'AntennaAngl2 ',self.AntennaAngl2
407 407 print 'RecPhaseCalibr0 ',self.RecPhaseCalibr0
408 408 print 'RecPhaseCalibr1 ',self.RecPhaseCalibr1
409 409 print 'RecPhaseCalibr2 ',self.RecPhaseCalibr2
410 410 print 'RecAmpCalibr0 ',self.RecAmpCalibr0
411 411 print 'RecAmpCalibr1 ',self.RecAmpCalibr1
412 412 print 'RecAmpCalibr2 ',self.RecAmpCalibr2
413 413 print 'ReceiverGaindB0 ',self.ReceiverGaindB0
414 414 print 'ReceiverGaindB1 ',self.ReceiverGaindB1
415 415 print 'ReceiverGaindB2 ',self.ReceiverGaindB2
416 416 print '=============================================='
417 417
418 418 if OffRHeader > endFp:
419 419 sys.stderr.write("Warning %s: Size value read from System Header is lower than it has to be\n" %fp)
420 420 return 0
421 421
422 422 if OffRHeader < endFp:
423 423 sys.stderr.write("Warning %s: Size value read from System Header size is greater than it has to be\n" %fp)
424 424 return 0
425 425
426 426 return 1
427 427
428 428
429 429 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODataReader):
430 430
431 431 path = None
432 432 startDate = None
433 433 endDate = None
434 434 startTime = None
435 435 endTime = None
436 436 walk = None
437 437 isConfig = False
438 438
439 439
440 440 fileList= None
441 441
442 442 #metadata
443 443 TimeZone= None
444 444 Interval= None
445 445 heightList= None
446 446
447 447 #data
448 448 data= None
449 449 utctime= None
450 450
451 451
452 452
453 453 def __init__(self, **kwargs):
454 454
455 455 #Eliminar de la base la herencia
456 456 ProcessingUnit.__init__(self, **kwargs)
457 457
458 458 #self.isConfig = False
459 459
460 460 #self.pts2read_SelfSpectra = 0
461 461 #self.pts2read_CrossSpectra = 0
462 462 #self.pts2read_DCchannels = 0
463 463 #self.datablock = None
464 464 self.utc = None
465 465 self.ext = ".fdt"
466 466 self.optchar = "P"
467 467 self.fpFile=None
468 468 self.fp = None
469 469 self.BlockCounter=0
470 470 self.dtype = None
471 471 self.fileSizeByHeader = None
472 472 self.filenameList = []
473 473 self.fileSelector = 0
474 474 self.Off2StartNxtRec=0
475 475 self.RecCounter=0
476 476 self.flagNoMoreFiles = 0
477 477 self.data_spc=None
478 478 self.data_cspc=None
479 479 self.data_output=None
480 480 self.path = None
481 481 self.OffsetStartHeader=0
482 482 self.Off2StartData=0
483 483 self.ipp = 0
484 484 self.nFDTdataRecors=0
485 485 self.blocksize = 0
486 486 self.dataOut = Spectra()
487 487 self.profileIndex = 1 #Always
488 488 self.dataOut.flagNoData=False
489 489 self.dataOut.nRdPairs = 0
490 490 self.dataOut.pairsList = []
491 491 self.dataOut.data_spc=None
492 492 self.dataOut.noise=[]
493 493 self.dataOut.velocityX=[]
494 494 self.dataOut.velocityY=[]
495 495 self.dataOut.velocityV=[]
496 496
497 497
498 498
499 499 def Files2Read(self, fp):
500 500 '''
501 501 Function that indicates the number of .fdt files that exist in the folder to be read.
502 502 It also creates an organized list with the names of the files to read.
503 503 '''
504 504 #self.__checkPath()
505 505
506 506 ListaData=os.listdir(fp) #Gets the list of files within the fp address
507 507 ListaData=sorted(ListaData) #Sort the list of files from least to largest by names
508 508 nFiles=0 #File Counter
509 509 FileList=[] #A list is created that will contain the .fdt files
510 510 for IndexFile in ListaData :
511 511 if '.fdt' in IndexFile:
512 512 FileList.append(IndexFile)
513 513 nFiles+=1
514 514
515 515 #print 'Files2Read'
516 516 #print 'Existen '+str(nFiles)+' archivos .fdt'
517 517
518 518 self.filenameList=FileList #List of files from least to largest by names
519 519
520 520
521 521 def run(self, **kwargs):
522 522 '''
523 523 This method will be the one that will initiate the data entry, will be called constantly.
524 524 You should first verify that your Setup () is set up and then continue to acquire
525 525 the data to be processed with getData ().
526 526 '''
527 527 if not self.isConfig:
528 528 self.setup(**kwargs)
529 529 self.isConfig = True
530 530
531 531 self.getData()
532 532 #print 'running'
533 533
534 534
535 535 def setup(self, path=None,
536 536 startDate=None,
537 537 endDate=None,
538 538 startTime=None,
539 539 endTime=None,
540 540 walk=True,
541 541 timezone='utc',
542 542 code = None,
543 543 online=False,
544 544 ReadMode=None,
545 545 **kwargs):
546 546
547 547 self.isConfig = True
548 548
549 549 self.path=path
550 550 self.startDate=startDate
551 551 self.endDate=endDate
552 552 self.startTime=startTime
553 553 self.endTime=endTime
554 554 self.walk=walk
555 555 self.ReadMode=int(ReadMode)
556 556
557 557 pass
558 558
559 559
560 560 def getData(self):
561 561 '''
562 562 Before starting this function, you should check that there is still an unread file,
563 563 If there are still blocks to read or if the data block is empty.
564 564
565 565 You should call the file "read".
566 566
567 567 '''
568 568
569 569 if self.flagNoMoreFiles:
570 570 self.dataOut.flagNoData = True
571 571 #print 'NoData se vuelve true'
572 572 return 0
573 573
574 574 self.fp=self.path
575 575 self.Files2Read(self.fp)
576 576 self.readFile(self.fp)
577 577 self.dataOut.data_spc = self.data_spc
578 578 self.dataOut.data_cspc =self.data_cspc
579 579 self.dataOut.data_output=self.data_output
580 580
581 581 #print 'self.dataOut.data_output', shape(self.dataOut.data_output)
582 582
583 583 #self.removeDC()
584 584 return self.dataOut.data_spc
585 585
586 586
587 587 def readFile(self,fp):
588 588 '''
589 589 You must indicate if you are reading in Online or Offline mode and load the
590 590 The parameters for this file reading mode.
591 591
592 592 Then you must do 2 actions:
593 593
594 594 1. Get the BLTR FileHeader.
595 595 2. Start reading the first block.
596 596 '''
597 597
598 598 #The address of the folder is generated the name of the .fdt file that will be read
599 599 #print "File: ",self.fileSelector+1
600 600
601 601 if self.fileSelector < len(self.filenameList):
602 602
603 603 self.fpFile=str(fp)+'/'+str(self.filenameList[self.fileSelector])
604 604 #print self.fpFile
605 605 fheader = FileHeaderBLTR()
606 606 fheader.FHread(self.fpFile) #Bltr FileHeader Reading
607 607 self.nFDTdataRecors=fheader.nFDTdataRecors
608 608
609 609 self.readBlock() #Block reading
610 610 else:
611 611 #print 'readFile FlagNoData becomes true'
612 612 self.flagNoMoreFiles=True
613 613 self.dataOut.flagNoData = True
614 614 return 0
615 615
616 616 def getVelRange(self, extrapoints=0):
617 617 Lambda= SPEED_OF_LIGHT/50000000
618 618 PRF = self.dataOut.PRF#1./(self.dataOut.ippSeconds * self.dataOut.nCohInt)
619 619 Vmax=-Lambda/(4.*(1./PRF)*self.dataOut.nCohInt*2.)
620 620 deltafreq = PRF / (self.nProfiles)
621 621 deltavel = (Vmax*2) / (self.nProfiles)
622 622 freqrange = deltafreq*(numpy.arange(self.nProfiles)-self.nProfiles/2.) - deltafreq/2
623 623 velrange = deltavel*(numpy.arange(self.nProfiles)-self.nProfiles/2.)
624 624 return velrange
625 625
626 626 def readBlock(self):
627 627 '''
628 628 It should be checked if the block has data, if it is not passed to the next file.
629 629
630 630 Then the following is done:
631 631
632 632 1. Read the RecordHeader
633 633 2. Fill the buffer with the current block number.
634 634
635 635 '''
636 636
637 637 if self.BlockCounter < self.nFDTdataRecors-1:
638 638 #print self.nFDTdataRecors, 'CONDICION'
639 639 if self.ReadMode==1:
640 640 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
641 641 elif self.ReadMode==0:
642 642 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
643 643
644 644 rheader.RHread(self.fpFile) #Bltr FileHeader Reading
645 645
646 646 self.OffsetStartHeader=rheader.OffsetStartHeader
647 647 self.RecCounter=rheader.RecCounter
648 648 self.Off2StartNxtRec=rheader.Off2StartNxtRec
649 649 self.Off2StartData=rheader.Off2StartData
650 650 self.nProfiles=rheader.nProfiles
651 651 self.nChannels=rheader.nChannels
652 652 self.nHeights=rheader.nHeights
653 653 self.frequency=rheader.TransmitFrec
654 654 self.DualModeIndex=rheader.DualModeIndex
655 655
656 656 self.pairsList =[(0,1),(0,2),(1,2)]
657 657 self.dataOut.pairsList = self.pairsList
658 658
659 659 self.nRdPairs=len(self.dataOut.pairsList)
660 660 self.dataOut.nRdPairs = self.nRdPairs
661 661
662 662 self.__firstHeigth=rheader.StartRangeSamp
663 663 self.__deltaHeigth=rheader.SampResolution
664 664 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
665 665 self.dataOut.channelList = range(self.nChannels)
666 666 self.dataOut.nProfiles=rheader.nProfiles
667 667 self.dataOut.nIncohInt=rheader.nIncohInt
668 668 self.dataOut.nCohInt=rheader.nCohInt
669 669 self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
670 670 self.dataOut.PRF=rheader.PRFhz
671 671 self.dataOut.nFFTPoints=rheader.nProfiles
672 672 self.dataOut.utctime=rheader.nUtime
673 673 self.dataOut.timeZone=0
674 674 self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt
675 675 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
676 676
677 677 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
678 678 #print 'self.data_output', shape(self.data_output)
679 679 self.dataOut.velocityX=[]
680 680 self.dataOut.velocityY=[]
681 681 self.dataOut.velocityV=[]
682 682
683 683 '''Block Reading, the Block Data is received and Reshape is used to give it
684 684 shape.
685 685 '''
686 686
687 687 #Procedure to take the pointer to where the date block starts
688 688 startDATA = open(self.fpFile,"rb")
689 689 OffDATA= self.OffsetStartHeader + self.RecCounter*self.Off2StartNxtRec+self.Off2StartData
690 690 startDATA.seek(OffDATA, os.SEEK_SET)
691 691
692 692 def moving_average(x, N=2):
693 693 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
694 694
695 695 def gaus(xSamples,a,x0,sigma):
696 696 return a*exp(-(xSamples-x0)**2/(2*sigma**2))
697 697
698 698 def Find(x,value):
699 699 for index in range(len(x)):
700 700 if x[index]==value:
701 701 return index
702 702
703 703 def pol2cart(rho, phi):
704 704 x = rho * numpy.cos(phi)
705 705 y = rho * numpy.sin(phi)
706 706 return(x, y)
707 707
708 708
709 709
710 710
711 711 if self.DualModeIndex==self.ReadMode:
712 712
713 713 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
714 714 #
715 715 # if len(self.data_fft) is not 101376:
716 716 #
717 717 # self.data_fft = numpy.empty(101376)
718 718
719 719 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
720 720
721 721
722 722
723 723
724 724
725 725 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
726 726
727 727 self.data_block = numpy.transpose(self.data_block, (1,2,0))
728 728
729 729 copy = self.data_block.copy()
730 730 spc = copy * numpy.conjugate(copy)
731 731
732 732 self.data_spc = numpy.absolute(spc) # valor absoluto o magnitud
733 733
734 734 factor = self.dataOut.normFactor
735 735
736 736
737 737 z = self.data_spc.copy()#/factor
738 738 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
739 739 #zdB = 10*numpy.log10(z)
740 740
741 741
742 742 self.dataOut.data_spc=self.data_spc
743 743
744 744 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
745 745 #noisedB = 10*numpy.log10(self.noise)
746 746
747 747
748 748 ySamples=numpy.ones([3,self.nProfiles])
749 749 phase=numpy.ones([3,self.nProfiles])
750 750 CSPCSamples=numpy.ones([3,self.nProfiles],dtype=numpy.complex_)
751 751 coherence=numpy.ones([3,self.nProfiles])
752 752 PhaseSlope=numpy.ones(3)
753 753 PhaseInter=numpy.ones(3)
754 754
755 755 '''****** Getting CrossSpectra ******'''
756 756 cspc=self.data_block.copy()
757 757 self.data_cspc=self.data_block.copy()
758 758
759 759 xFrec=self.getVelRange(1)
760 760 VelRange=self.getVelRange(1)
761 761 self.dataOut.VelRange=VelRange
762 762 #print ' '
763 763 #print ' '
764 764 #print 'xFrec',xFrec
765 765 #print ' '
766 766 #print ' '
767 767 #Height=35
768 768
769 769 for i in range(self.nRdPairs):
770 770
771 771 chan_index0 = self.dataOut.pairsList[i][0]
772 772 chan_index1 = self.dataOut.pairsList[i][1]
773 773
774 774 self.data_cspc[i,:,:]=cspc[chan_index0,:,:] * numpy.conjugate(cspc[chan_index1,:,:])
775 775
776 776
777 777 '''Getting Eij and Nij'''
778 778 (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
779 779 (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
780 780 (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
781 781
782 782 E01=AntennaX0-AntennaX1
783 783 N01=AntennaY0-AntennaY1
784 784
785 785 E02=AntennaX0-AntennaX2
786 786 N02=AntennaY0-AntennaY2
787 787
788 788 E12=AntennaX1-AntennaX2
789 789 N12=AntennaY1-AntennaY2
790 790
791 791 self.ChanDist= numpy.array([[E01, N01],[E02,N02],[E12,N12]])
792 792
793 793 self.dataOut.ChanDist = self.ChanDist
794 794
795 795
796
797
798 self.BlockCounter+=2
799
800 else:
801 self.fileSelector+=1
802 self.BlockCounter=0
803 print "Next File" No newline at end of file
@@ -1,3809 +1,3811
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import scipy
5 5 import re
6 6 import datetime
7 7 import copy
8 8 import sys
9 9 import importlib
10 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 13
14 14
15 15
16 16 import time
17 17 #from sklearn.cluster import KMeans
18 18
19 19 import matplotlib.pyplot as plt
20 20
21 21 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
22 22 from jroproc_base import ProcessingUnit, Operation
23 23 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
24 24 from scipy import asarray as ar,exp
25 25 from scipy.optimize import curve_fit
26 26
27 27 import warnings
28 28 from numpy import NaN
29 29 from scipy.optimize.optimize import OptimizeWarning
30 30 warnings.filterwarnings('ignore')
31 31
32 32
33 33 SPEED_OF_LIGHT = 299792458
34 34
35 35
36 36 '''solving pickling issue'''
37 37
38 38 def _pickle_method(method):
39 39 func_name = method.im_func.__name__
40 40 obj = method.im_self
41 41 cls = method.im_class
42 42 return _unpickle_method, (func_name, obj, cls)
43 43
44 44 def _unpickle_method(func_name, obj, cls):
45 45 for cls in cls.mro():
46 46 try:
47 47 func = cls.__dict__[func_name]
48 48 except KeyError:
49 49 pass
50 50 else:
51 51 break
52 52 return func.__get__(obj, cls)
53 53
54 54 class ParametersProc(ProcessingUnit):
55 55
56 56 nSeconds = None
57 57
58 58 def __init__(self):
59 59 ProcessingUnit.__init__(self)
60 60
61 61 # self.objectDict = {}
62 62 self.buffer = None
63 63 self.firstdatatime = None
64 64 self.profIndex = 0
65 65 self.dataOut = Parameters()
66 66
67 67 def __updateObjFromInput(self):
68 68
69 69 self.dataOut.inputUnit = self.dataIn.type
70 70
71 71 self.dataOut.timeZone = self.dataIn.timeZone
72 72 self.dataOut.dstFlag = self.dataIn.dstFlag
73 73 self.dataOut.errorCount = self.dataIn.errorCount
74 74 self.dataOut.useLocalTime = self.dataIn.useLocalTime
75 75
76 76 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
77 77 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
78 78 self.dataOut.channelList = self.dataIn.channelList
79 79 self.dataOut.heightList = self.dataIn.heightList
80 80 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
81 81 # self.dataOut.nHeights = self.dataIn.nHeights
82 82 # self.dataOut.nChannels = self.dataIn.nChannels
83 83 self.dataOut.nBaud = self.dataIn.nBaud
84 84 self.dataOut.nCode = self.dataIn.nCode
85 85 self.dataOut.code = self.dataIn.code
86 86 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
87 87 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
88 88 # self.dataOut.utctime = self.firstdatatime
89 89 self.dataOut.utctime = self.dataIn.utctime
90 90 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
91 91 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
92 92 self.dataOut.nCohInt = self.dataIn.nCohInt
93 93 # self.dataOut.nIncohInt = 1
94 94 self.dataOut.ippSeconds = self.dataIn.ippSeconds
95 95 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
96 96 self.dataOut.timeInterval1 = self.dataIn.timeInterval
97 97 self.dataOut.heightList = self.dataIn.getHeiRange()
98 98 self.dataOut.frequency = self.dataIn.frequency
99 99 self.dataOut.noise = self.dataIn.noise
100 100
101 101
102 102
103 103 def run(self):
104 104
105 105 #---------------------- Voltage Data ---------------------------
106 106
107 107 if self.dataIn.type == "Voltage":
108 108
109 109 self.__updateObjFromInput()
110 110 self.dataOut.data_pre = self.dataIn.data.copy()
111 111 self.dataOut.flagNoData = False
112 112 self.dataOut.utctimeInit = self.dataIn.utctime
113 113 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
114 114 return
115 115
116 116 #---------------------- Spectra Data ---------------------------
117 117
118 118 if self.dataIn.type == "Spectra":
119 119
120 120 self.dataOut.data_pre = (self.dataIn.data_spc , self.dataIn.data_cspc)
121 121 print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
122 122 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
123 123 self.dataOut.spc_noise = self.dataIn.getNoise()
124 124 self.dataOut.spc_range = numpy.asanyarray((self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) ))
125 125
126 126 self.dataOut.normFactor = self.dataIn.normFactor
127 127 #self.dataOut.outputInterval = self.dataIn.outputInterval
128 128 self.dataOut.groupList = self.dataIn.pairsList
129 129 self.dataOut.flagNoData = False
130 130 #print 'datain chandist ',self.dataIn.ChanDist
131 131 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
132 132 self.dataOut.ChanDist = self.dataIn.ChanDist
133 133 else: self.dataOut.ChanDist = None
134 134
135 135 print 'datain chandist ',self.dataOut.ChanDist
136 136
137 137 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
138 138 # self.dataOut.VelRange = self.dataIn.VelRange
139 139 #else: self.dataOut.VelRange = None
140 140
141 141 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
142 142 self.dataOut.RadarConst = self.dataIn.RadarConst
143 143
144 144 if hasattr(self.dataIn, 'NPW'): #NPW
145 145 self.dataOut.NPW = self.dataIn.NPW
146 146
147 147 if hasattr(self.dataIn, 'COFA'): #COFA
148 148 self.dataOut.COFA = self.dataIn.COFA
149 149
150 150
151 151
152 152 #---------------------- Correlation Data ---------------------------
153 153
154 154 if self.dataIn.type == "Correlation":
155 155 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
156 156
157 157 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
158 158 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
159 159 self.dataOut.groupList = (acf_pairs, ccf_pairs)
160 160
161 161 self.dataOut.abscissaList = self.dataIn.lagRange
162 162 self.dataOut.noise = self.dataIn.noise
163 163 self.dataOut.data_SNR = self.dataIn.SNR
164 164 self.dataOut.flagNoData = False
165 165 self.dataOut.nAvg = self.dataIn.nAvg
166 166
167 167 #---------------------- Parameters Data ---------------------------
168 168
169 169 if self.dataIn.type == "Parameters":
170 170 self.dataOut.copy(self.dataIn)
171 171 self.dataOut.flagNoData = False
172 172
173 173 return True
174 174
175 175 self.__updateObjFromInput()
176 176 self.dataOut.utctimeInit = self.dataIn.utctime
177 177 self.dataOut.paramInterval = self.dataIn.timeInterval
178 178
179 179 return
180 180
181 181
182 182 def target(tups):
183 183
184 184 obj, args = tups
185 185
186 186 return obj.FitGau(args)
187 187
188 188
189 189 class SpectralFilters(Operation):
190 190
191 191 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
192 192
193 193 LimitR : It is the limit in m/s of Rainfall
194 194 LimitW : It is the limit in m/s for Winds
195 195
196 196 Input:
197 197
198 198 self.dataOut.data_pre : SPC and CSPC
199 199 self.dataOut.spc_range : To select wind and rainfall velocities
200 200
201 201 Affected:
202 202
203 203 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
204 204 self.dataOut.spcparam_range : Used in SpcParamPlot
205 205 self.dataOut.SPCparam : Used in PrecipitationProc
206 206
207 207
208 208 '''
209 209
210 210 def __init__(self, **kwargs):
211 211 Operation.__init__(self, **kwargs)
212 212 self.i=0
213 213
214 214 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
215 215
216 216
217 217 #Limite de vientos
218 218 LimitR = PositiveLimit
219 219 LimitN = NegativeLimit
220 220
221 221 self.spc = dataOut.data_pre[0].copy()
222 222 self.cspc = dataOut.data_pre[1].copy()
223 223
224 224 self.Num_Hei = self.spc.shape[2]
225 225 self.Num_Bin = self.spc.shape[1]
226 226 self.Num_Chn = self.spc.shape[0]
227 227
228 228 VelRange = dataOut.spc_range[2]
229 229 TimeRange = dataOut.spc_range[1]
230 230 FrecRange = dataOut.spc_range[0]
231 231
232 232 Vmax= 2*numpy.max(dataOut.spc_range[2])
233 233 Tmax= 2*numpy.max(dataOut.spc_range[1])
234 234 Fmax= 2*numpy.max(dataOut.spc_range[0])
235 235
236 236 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
237 237 Breaker1R=numpy.where(VelRange == Breaker1R)
238 238
239 239 Delta = self.Num_Bin/2 - Breaker1R[0]
240 240
241 241
242 242 '''Reacomodando SPCrange'''
243 243
244 244 VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0)
245 245
246 246 VelRange[-(self.Num_Bin/2):]+= Vmax
247 247
248 248 FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0)
249 249
250 250 FrecRange[-(self.Num_Bin/2):]+= Fmax
251 251
252 252 TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0)
253 253
254 254 TimeRange[-(self.Num_Bin/2):]+= Tmax
255 255
256 256 ''' ------------------ '''
257 257
258 258 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
259 259 Breaker2R=numpy.where(VelRange == Breaker2R)
260 260
261 261
262 262 SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1)
263 263
264 264 SPCcut = SPCroll.copy()
265 265 for i in range(self.Num_Chn):
266 266
267 267 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
268 268 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
269 269
270 270 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
271 271 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
272 272
273 273 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
274 274 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
275 275
276 276 SPC_ch1 = SPCroll
277 277
278 278 SPC_ch2 = SPCcut
279 279
280 280 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
281 281 dataOut.SPCparam = numpy.asarray(SPCparam)
282 282
283 283
284 284 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
285 285
286 286 dataOut.spcparam_range[2]=VelRange
287 287 dataOut.spcparam_range[1]=TimeRange
288 288 dataOut.spcparam_range[0]=FrecRange
289 289
290 290
291 291 class GaussianFit(Operation):
292 292
293 293 '''
294 294 Function that fit of one and two generalized gaussians (gg) based
295 295 on the PSD shape across an "power band" identified from a cumsum of
296 296 the measured spectrum - noise.
297 297
298 298 Input:
299 299 self.dataOut.data_pre : SelfSpectra
300 300
301 301 Output:
302 302 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
303 303
304 304 '''
305 305 def __init__(self, **kwargs):
306 306 Operation.__init__(self, **kwargs)
307 307 self.i=0
308 308
309 309
310 310 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
311 311 """This routine will find a couple of generalized Gaussians to a power spectrum
312 312 input: spc
313 313 output:
314 314 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
315 315 """
316 316
317 317 self.spc = dataOut.data_pre[0].copy()
318 318
319 319
320 320 self.Num_Hei = self.spc.shape[2]
321 321 self.Num_Bin = self.spc.shape[1]
322 322 self.Num_Chn = self.spc.shape[0]
323 323 Vrange = dataOut.abscissaList
324 324
325 325 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
326 326 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
327 327 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
328 328 SPC_ch1[:] = numpy.NaN
329 329 SPC_ch2[:] = numpy.NaN
330 330
331 331
332 332 start_time = time.time()
333 333
334 334 noise_ = dataOut.spc_noise[0].copy()
335 335
336 336
337 337 pool = Pool(processes=self.Num_Chn)
338 338 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
339 339 objs = [self for __ in range(self.Num_Chn)]
340 340 attrs = zip(objs, args)
341 341 gauSPC = pool.map(target, attrs)
342 342 dataOut.SPCparam = numpy.asarray(SPCparam)
343 343
344 344 # re-normalizing spc and noise
345 345 # This part differs from gg1
346 346
347 347 ''' Parameters:
348 348 1. Amplitude
349 349 2. Shift
350 350 3. Width
351 351 4. Power
352 352 '''
353 353
354 354
355 355 ###############################################################################
356 356 def FitGau(self, X):
357 357
358 358 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
359 359
360 360 SPCparam = []
361 361 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
362 362 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
363 363 SPC_ch1[:] = 0#numpy.NaN
364 364 SPC_ch2[:] = 0#numpy.NaN
365 365
366 366
367 367
368 368 for ht in range(self.Num_Hei):
369 369
370 370
371 371 spc = numpy.asarray(self.spc)[ch,:,ht]
372 372
373 373 #############################################
374 374 # normalizing spc and noise
375 375 # This part differs from gg1
376 376 spc_norm_max = max(spc)
377 377 #spc = spc / spc_norm_max
378 378 pnoise = pnoise #/ spc_norm_max
379 379 #############################################
380 380
381 381 fatspectra=1.0
382 382
383 383 wnoise = noise_ #/ spc_norm_max
384 384 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
385 385 #if wnoise>1.1*pnoise: # to be tested later
386 386 # wnoise=pnoise
387 387 noisebl=wnoise*0.9;
388 388 noisebh=wnoise*1.1
389 389 spc=spc-wnoise
390 390
391 391 minx=numpy.argmin(spc)
392 392 #spcs=spc.copy()
393 393 spcs=numpy.roll(spc,-minx)
394 394 cum=numpy.cumsum(spcs)
395 395 tot_noise=wnoise * self.Num_Bin #64;
396 396
397 397 snr = sum(spcs)/tot_noise
398 398 snrdB=10.*numpy.log10(snr)
399 399
400 400 if snrdB < SNRlimit :
401 401 snr = numpy.NaN
402 402 SPC_ch1[:,ht] = 0#numpy.NaN
403 403 SPC_ch1[:,ht] = 0#numpy.NaN
404 404 SPCparam = (SPC_ch1,SPC_ch2)
405 405 continue
406 406
407 407
408 408 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
409 409 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
410 410
411 411 cummax=max(cum);
412 412 epsi=0.08*fatspectra # cumsum to narrow down the energy region
413 413 cumlo=cummax*epsi;
414 414 cumhi=cummax*(1-epsi)
415 415 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
416 416
417 417
418 418 if len(powerindex) < 1:# case for powerindex 0
419 419 continue
420 420 powerlo=powerindex[0]
421 421 powerhi=powerindex[-1]
422 422 powerwidth=powerhi-powerlo
423 423
424 424 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
425 425 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
426 426 midpeak=(firstpeak+secondpeak)/2.
427 427 firstamp=spcs[int(firstpeak)]
428 428 secondamp=spcs[int(secondpeak)]
429 429 midamp=spcs[int(midpeak)]
430 430
431 431 x=numpy.arange( self.Num_Bin )
432 432 y_data=spc+wnoise
433 433
434 434 ''' single Gaussian '''
435 435 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
436 436 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
437 437 power0=2.
438 438 amplitude0=midamp
439 439 state0=[shift0,width0,amplitude0,power0,wnoise]
440 440 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
441 441 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
442 442
443 443 chiSq1=lsq1[1];
444 444
445 445
446 446 if fatspectra<1.0 and powerwidth<4:
447 447 choice=0
448 448 Amplitude0=lsq1[0][2]
449 449 shift0=lsq1[0][0]
450 450 width0=lsq1[0][1]
451 451 p0=lsq1[0][3]
452 452 Amplitude1=0.
453 453 shift1=0.
454 454 width1=0.
455 455 p1=0.
456 456 noise=lsq1[0][4]
457 457 #return (numpy.array([shift0,width0,Amplitude0,p0]),
458 458 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
459 459
460 460 ''' two gaussians '''
461 461 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
462 462 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
463 463 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
464 464 width0=powerwidth/6.;
465 465 width1=width0
466 466 power0=2.;
467 467 power1=power0
468 468 amplitude0=firstamp;
469 469 amplitude1=secondamp
470 470 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
471 471 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
472 472 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
473 473 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
474 474
475 475 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
476 476
477 477
478 478 chiSq2=lsq2[1];
479 479
480 480
481 481
482 482 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
483 483
484 484 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
485 485 if oneG:
486 486 choice=0
487 487 else:
488 488 w1=lsq2[0][1]; w2=lsq2[0][5]
489 489 a1=lsq2[0][2]; a2=lsq2[0][6]
490 490 p1=lsq2[0][3]; p2=lsq2[0][7]
491 491 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
492 492 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
493 493 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
494 494
495 495 if gp1>gp2:
496 496 if a1>0.7*a2:
497 497 choice=1
498 498 else:
499 499 choice=2
500 500 elif gp2>gp1:
501 501 if a2>0.7*a1:
502 502 choice=2
503 503 else:
504 504 choice=1
505 505 else:
506 506 choice=numpy.argmax([a1,a2])+1
507 507 #else:
508 508 #choice=argmin([std2a,std2b])+1
509 509
510 510 else: # with low SNR go to the most energetic peak
511 511 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
512 512
513 513
514 514 shift0=lsq2[0][0];
515 515 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
516 516 shift1=lsq2[0][4];
517 517 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
518 518
519 519 max_vel = 1.0
520 520
521 521 #first peak will be 0, second peak will be 1
522 522 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
523 523 shift0=lsq2[0][0]
524 524 width0=lsq2[0][1]
525 525 Amplitude0=lsq2[0][2]
526 526 p0=lsq2[0][3]
527 527
528 528 shift1=lsq2[0][4]
529 529 width1=lsq2[0][5]
530 530 Amplitude1=lsq2[0][6]
531 531 p1=lsq2[0][7]
532 532 noise=lsq2[0][8]
533 533 else:
534 534 shift1=lsq2[0][0]
535 535 width1=lsq2[0][1]
536 536 Amplitude1=lsq2[0][2]
537 537 p1=lsq2[0][3]
538 538
539 539 shift0=lsq2[0][4]
540 540 width0=lsq2[0][5]
541 541 Amplitude0=lsq2[0][6]
542 542 p0=lsq2[0][7]
543 543 noise=lsq2[0][8]
544 544
545 545 if Amplitude0<0.05: # in case the peak is noise
546 546 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
547 547 if Amplitude1<0.05:
548 548 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
549 549
550 550
551 551 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
552 552 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
553 553 SPCparam = (SPC_ch1,SPC_ch2)
554 554
555 555
556 556 return GauSPC
557 557
558 558 def y_model1(self,x,state):
559 559 shift0,width0,amplitude0,power0,noise=state
560 560 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
561 561
562 562 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
563 563
564 564 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
565 565 return model0+model0u+model0d+noise
566 566
567 567 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
568 568 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
569 569 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
570 570
571 571 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
572 572
573 573 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
574 574 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
575 575
576 576 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
577 577
578 578 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
579 579 return model0+model0u+model0d+model1+model1u+model1d+noise
580 580
581 581 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
582 582
583 583 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
584 584
585 585 def misfit2(self,state,y_data,x,num_intg):
586 586 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
587 587
588 588
589 589
590 590 class PrecipitationProc(Operation):
591 591
592 592 '''
593 593 Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
594 594
595 595 Input:
596 596 self.dataOut.data_pre : SelfSpectra
597 597
598 598 Output:
599 599
600 600 self.dataOut.data_output : Reflectivity factor, rainfall Rate
601 601
602 602
603 603 Parameters affected:
604 604 '''
605 605
606 606 def __init__(self, **kwargs):
607 607 Operation.__init__(self, **kwargs)
608 608 self.i=0
609 609
610 610
611 611 def gaus(self,xSamples,Amp,Mu,Sigma):
612 612 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
613 613
614 614
615 615
616 616 def Moments(self, ySamples, xSamples):
617 617 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
618 618 yNorm = ySamples / Pot
619 619
620 620 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
621 621 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
622 622 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
623 623
624 624 return numpy.array([Pot, Vr, Desv])
625 625
626 626 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
627 627 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
628 628
629 629
630 630 Velrange = dataOut.spcparam_range[2]
631 631 FrecRange = dataOut.spcparam_range[0]
632 632
633 633 dV= Velrange[1]-Velrange[0]
634 634 dF= FrecRange[1]-FrecRange[0]
635 635
636 636 if radar == "MIRA35C" :
637 637
638 638 self.spc = dataOut.data_pre[0].copy()
639 639 self.Num_Hei = self.spc.shape[2]
640 640 self.Num_Bin = self.spc.shape[1]
641 641 self.Num_Chn = self.spc.shape[0]
642 642 Ze = self.dBZeMODE2(dataOut)
643 643
644 644 else:
645 645
646 646 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
647 647
648 648 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
649 649
650 650 self.spc[:,:,0:7]= numpy.NaN
651 651
652 652 """##########################################"""
653 653
654 654 self.Num_Hei = self.spc.shape[2]
655 655 self.Num_Bin = self.spc.shape[1]
656 656 self.Num_Chn = self.spc.shape[0]
657 657
658 658 ''' Se obtiene la constante del RADAR '''
659 659
660 660 self.Pt = Pt
661 661 self.Gt = Gt
662 662 self.Gr = Gr
663 663 self.Lambda = Lambda
664 664 self.aL = aL
665 665 self.tauW = tauW
666 666 self.ThetaT = ThetaT
667 667 self.ThetaR = ThetaR
668 668
669 669 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
670 670 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
671 671 RadarConstant = 5e-26 * Numerator / Denominator #
672 672
673 673 ''' ============================= '''
674 674
675 675 self.spc[0] = (self.spc[0]-dataOut.noise[0])
676 676 self.spc[1] = (self.spc[1]-dataOut.noise[1])
677 677 self.spc[2] = (self.spc[2]-dataOut.noise[2])
678 678
679 679 self.spc[ numpy.where(self.spc < 0)] = 0
680 680
681 681 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
682 682 SPCmean[ numpy.where(SPCmean < 0)] = 0
683 683
684 684 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
685 685 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
686 686 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
687 687
688 688 Pr = SPCmean[:,:]
689 689
690 690 VelMeteoro = numpy.mean(SPCmean,axis=0)
691 691
692 692 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
693 693 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
694 694 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
695 695 V_mean = numpy.zeros(self.Num_Hei)
696 696 del_V = numpy.zeros(self.Num_Hei)
697 697 Z = numpy.zeros(self.Num_Hei)
698 698 Ze = numpy.zeros(self.Num_Hei)
699 699 RR = numpy.zeros(self.Num_Hei)
700 700
701 701 Range = dataOut.heightList*1000.
702 702
703 703 for R in range(self.Num_Hei):
704 704
705 705 h = Range[R] + Altitude #Range from ground to radar pulse altitude
706 706 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
707 707
708 708 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
709 709
710 710 '''NOTA: ETA(n) dn = ETA(f) df
711 711
712 712 dn = 1 Diferencial de muestreo
713 713 df = ETA(n) / ETA(f)
714 714
715 715 '''
716 716
717 717 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
718 718
719 719 ETAv[:,R]=ETAn[:,R]/dV
720 720
721 721 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
722 722
723 723 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
724 724
725 725 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
726 726
727 727 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
728 728
729 729 try:
730 730 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
731 731 except:
732 732 popt01=numpy.zeros(3)
733 733 popt01[1]= DMoments[1]
734 734
735 735 if popt01[1]<0 or popt01[1]>20:
736 736 popt01[1]=numpy.NaN
737 737
738 738
739 739 V_mean[R]=popt01[1]
740 740
741 741 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
742 742
743 743 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
744 744
745 745 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
746 746
747 747
748 748
749 749 RR2 = (Z/200)**(1/1.6)
750 750 dBRR = 10*numpy.log10(RR)
751 751 dBRR2 = 10*numpy.log10(RR2)
752 752
753 753 dBZe = 10*numpy.log10(Ze)
754 754 dBZ = 10*numpy.log10(Z)
755 755
756 756 dataOut.data_output = RR[8]
757 757 dataOut.data_param = numpy.ones([3,self.Num_Hei])
758 758 dataOut.channelList = [0,1,2]
759 759
760 760 dataOut.data_param[0]=dBZ
761 761 dataOut.data_param[1]=V_mean
762 762 dataOut.data_param[2]=RR
763 763
764 764
765 765 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
766 766
767 767 NPW = dataOut.NPW
768 768 COFA = dataOut.COFA
769 769
770 770 SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
771 771 RadarConst = dataOut.RadarConst
772 772 #frequency = 34.85*10**9
773 773
774 774 ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
775 775 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
776 776
777 777 ETA = numpy.sum(SNR,1)
778 778
779 779 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
780 780
781 781 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
782 782
783 783 for r in range(self.Num_Hei):
784 784
785 785 Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
786 786 #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
787 787
788 788 return Ze
789 789
790 790 # def GetRadarConstant(self):
791 791 #
792 792 # """
793 793 # Constants:
794 794 #
795 795 # Pt: Transmission Power dB 5kW 5000
796 796 # Gt: Transmission Gain dB 24.7 dB 295.1209
797 797 # Gr: Reception Gain dB 18.5 dB 70.7945
798 798 # Lambda: Wavelenght m 0.6741 m 0.6741
799 799 # aL: Attenuation loses dB 4dB 2.5118
800 800 # tauW: Width of transmission pulse s 4us 4e-6
801 801 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
802 802 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
803 803 #
804 804 # """
805 805 #
806 806 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
807 807 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
808 808 # RadarConstant = Numerator / Denominator
809 809 #
810 810 # return RadarConstant
811 811
812 812
813 813
814 814 class FullSpectralAnalysis(Operation):
815 815
816 816 """
817 817 Function that implements Full Spectral Analisys technique.
818 818
819 819 Input:
820 820 self.dataOut.data_pre : SelfSpectra and CrossSPectra data
821 821 self.dataOut.groupList : Pairlist of channels
822 822 self.dataOut.ChanDist : Physical distance between receivers
823 823
824 824
825 825 Output:
826 826
827 827 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
828 828
829 829
830 830 Parameters affected: Winds, height range, SNR
831 831
832 832 """
833 833 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7):
834 834
835 835 self.indice=int(numpy.random.rand()*1000)
836 836
837 837 spc = dataOut.data_pre[0].copy()
838 838 cspc = dataOut.data_pre[1]
839 839
840 840 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
841 841
842 842 SNRspc = spc.copy()
843 843 SNRspc[:,:,0:7]= numpy.NaN
844 844
845 845 """##########################################"""
846 846
847 847
848 848 nChannel = spc.shape[0]
849 849 nProfiles = spc.shape[1]
850 850 nHeights = spc.shape[2]
851 851
852 852 pairsList = dataOut.groupList
853 853 if dataOut.ChanDist is not None :
854 854 ChanDist = dataOut.ChanDist
855 855 else:
856 856 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
857 857
858 858 FrecRange = dataOut.spc_range[0]
859 859
860 860 ySamples=numpy.ones([nChannel,nProfiles])
861 861 phase=numpy.ones([nChannel,nProfiles])
862 862 CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
863 863 coherence=numpy.ones([nChannel,nProfiles])
864 864 PhaseSlope=numpy.ones(nChannel)
865 865 PhaseInter=numpy.ones(nChannel)
866 866 data_SNR=numpy.zeros([nProfiles])
867 867
868 868 data = dataOut.data_pre
869 869 noise = dataOut.noise
870 870
871 871 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
872 872
873 873 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
874 874
875 875
876 876 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
877 877
878 878 velocityX=[]
879 879 velocityY=[]
880 880 velocityV=[]
881 881 PhaseLine=[]
882 882
883 883 dbSNR = 10*numpy.log10(dataOut.data_SNR)
884 884 dbSNR = numpy.average(dbSNR,0)
885 885
886 886 for Height in range(nHeights):
887 887
888 888 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit)
889 889 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
890 890
891 891 if abs(Vzon)<100. and abs(Vzon)> 0.:
892 892 velocityX=numpy.append(velocityX, Vzon)#Vmag
893 893
894 894 else:
895 895 velocityX=numpy.append(velocityX, numpy.NaN)
896 896
897 897 if abs(Vmer)<100. and abs(Vmer) > 0.:
898 898 velocityY=numpy.append(velocityY, -Vmer)#Vang
899 899
900 900 else:
901 901 velocityY=numpy.append(velocityY, numpy.NaN)
902 902
903 903 if dbSNR[Height] > SNRlimit:
904 904 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
905 905 else:
906 906 velocityV=numpy.append(velocityV, numpy.NaN)
907 907
908 908
909 909
910
910 '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR'''
911 911 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
912 912 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
913 913 data_output[2] = velocityV#FirstMoment
914 914
915 915 xFrec=FrecRange[0:spc.shape[1]]
916 916
917 917 dataOut.data_output=data_output
918 918
919 print 'velocidad zonal', data_output[0]
920
919 921 return
920 922
921 923
922 924 def moving_average(self,x, N=2):
923 925 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
924 926
925 927 def gaus(self,xSamples,Amp,Mu,Sigma):
926 928 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
927 929
928 930
929 931
930 932 def Moments(self, ySamples, xSamples):
931 933 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
932 934 yNorm = ySamples / Pot
933 935 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
934 936 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
935 937 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
936 938
937 939 return numpy.array([Pot, Vr, Desv])
938 940
939 941 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
940 942
941 943
942 944
943 945 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
944 946 phase=numpy.ones([spc.shape[0],spc.shape[1]])
945 947 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
946 948 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
947 949 PhaseSlope=numpy.zeros(spc.shape[0])
948 950 PhaseInter=numpy.ones(spc.shape[0])
949 951 xFrec=AbbsisaRange[0][0:spc.shape[1]]
950 952 xVel =AbbsisaRange[2][0:spc.shape[1]]
951 953 Vv=numpy.empty(spc.shape[2])*0
952 954 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
953 955
954 956 SPCmoments = self.Moments(SPCav[:,Height], xVel )
955 957 CSPCmoments = []
956 958 cspcNoise = numpy.empty(3)
957 959
958 960 '''Getting Eij and Nij'''
959 961
960 962 Xi01=ChanDist[0][0]
961 963 Eta01=ChanDist[0][1]
962 964
963 965 Xi02=ChanDist[1][0]
964 966 Eta02=ChanDist[1][1]
965 967
966 968 Xi12=ChanDist[2][0]
967 969 Eta12=ChanDist[2][1]
968 970
969 971 z = spc.copy()
970 972 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
971 973
972 974 for i in range(spc.shape[0]):
973 975
974 976 '''****** Line of Data SPC ******'''
975 977 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
976 978
977 979 '''****** SPC is normalized ******'''
978 980 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
979 981 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
980 982
981 983 xSamples = xFrec # Se toma el rango de frecuncias
982 984 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
983 985
984 986 for i in range(spc.shape[0]):
985 987
986 988 '''****** Line of Data CSPC ******'''
987 989 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
988 990 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
989 991 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
990 992
991 993 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
992 994 chan_index0 = pairsList[i][0]
993 995 chan_index1 = pairsList[i][1]
994 996
995 997 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
996 998 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
997 999
998 1000 CSPCSamples[i] = CSPCNorm
999 1001
1000 1002 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1001 1003
1002 1004 #coherence[i]= self.moving_average(coherence[i],N=1)
1003 1005
1004 1006 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1005 1007
1006 1008 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
1007 1009 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1008 1010 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1009 1011
1010 1012
1011 1013 popt=[1e-10,0,1e-10]
1012 1014 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1013 1015 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1014 1016
1015 1017 CSPCMask01 = numpy.abs(CSPCSamples[0])
1016 1018 CSPCMask02 = numpy.abs(CSPCSamples[1])
1017 1019 CSPCMask12 = numpy.abs(CSPCSamples[2])
1018 1020
1019 1021 mask01 = ~numpy.isnan(CSPCMask01)
1020 1022 mask02 = ~numpy.isnan(CSPCMask02)
1021 1023 mask12 = ~numpy.isnan(CSPCMask12)
1022 1024
1023 1025 #mask = ~numpy.isnan(CSPCMask01)
1024 1026 CSPCMask01 = CSPCMask01[mask01]
1025 1027 CSPCMask02 = CSPCMask02[mask02]
1026 1028 CSPCMask12 = CSPCMask12[mask12]
1027 1029 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1028 1030
1029 1031
1030 1032
1031 1033 '''***Fit Gauss CSPC01***'''
1032 1034 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1033 1035 try:
1034 1036 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1035 1037 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1036 1038 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1037 1039 FitGauss01 = self.gaus(xSamples,*popt01)
1038 1040 FitGauss02 = self.gaus(xSamples,*popt02)
1039 1041 FitGauss12 = self.gaus(xSamples,*popt12)
1040 1042 except:
1041 1043 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
1042 1044 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
1043 1045 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
1044 1046
1045 1047
1046 1048 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1047 1049
1048 1050 '''****** Getting fij width ******'''
1049 1051
1050 1052 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1051 1053
1052 1054 '''******* Getting fitting Gaussian *******'''
1053 1055 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1054 1056 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1055 1057
1056 1058 yMoments = self.Moments(yMean, xSamples)
1057 1059
1058 1060 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1059 1061 try:
1060 1062 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1061 1063 FitGauss=self.gaus(xSamples,*popt)
1062 1064
1063 1065 except :#RuntimeError:
1064 1066 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1065 1067
1066 1068
1067 1069 else:
1068 1070 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1069 1071
1070 1072
1071 1073
1072 1074 '''****** Getting Fij ******'''
1073 1075 Fijcspc = CSPCopt[:,2]/2*3
1074 1076
1075 1077
1076 1078 GaussCenter = popt[1] #xFrec[GCpos]
1077 1079 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1078 1080 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1079 1081 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1080 1082
1081 1083 #Punto e^-1 hubicado en la Gaussiana
1082 1084 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1083 1085 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1084 1086 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1085 1087
1086 1088 if xSamples[PointFij] > xSamples[PointGauCenter]:
1087 1089 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1088 1090
1089 1091 else:
1090 1092 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1091 1093
1092 1094
1093 1095 '''****** Taking frequency ranges from SPCs ******'''
1094 1096
1095 1097
1096 1098 #GaussCenter = popt[1] #Primer momento 01
1097 1099 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1098 1100 Range = numpy.empty(2)
1099 1101 Range[0] = GaussCenter - GauWidth
1100 1102 Range[1] = GaussCenter + GauWidth
1101 1103 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1102 1104 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1103 1105 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1104 1106
1105 1107 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1106 1108 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1107 1109
1108 1110 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1109 1111
1110 1112 FrecRange = xFrec[ Range[0] : Range[1] ]
1111 1113 VelRange = xVel[ Range[0] : Range[1] ]
1112 1114
1113 1115
1114 1116 '''****** Getting SCPC Slope ******'''
1115 1117
1116 1118 for i in range(spc.shape[0]):
1117 1119
1118 1120 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1119 1121 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1120 1122
1121 1123 '''***********************VelRange******************'''
1122 1124
1123 1125 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1124 1126
1125 1127 if len(FrecRange) == len(PhaseRange):
1126 1128 try:
1127 1129 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1128 1130 PhaseSlope[i]=slope
1129 1131 PhaseInter[i]=intercept
1130 1132 except:
1131 1133 PhaseSlope[i]=0
1132 1134 PhaseInter[i]=0
1133 1135 else:
1134 1136 PhaseSlope[i]=0
1135 1137 PhaseInter[i]=0
1136 1138 else:
1137 1139 PhaseSlope[i]=0
1138 1140 PhaseInter[i]=0
1139 1141
1140 1142
1141 1143 '''Getting constant C'''
1142 1144 cC=(Fij*numpy.pi)**2
1143 1145
1144 1146 '''****** Getting constants F and G ******'''
1145 1147 MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1146 1148 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1147 1149 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1148 1150 MijResults=numpy.array([MijResult0,MijResult1])
1149 1151 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1150 1152
1151 1153 '''****** Getting constants A, B and H ******'''
1152 1154 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1153 1155 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1154 1156 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1155 1157
1156 1158 WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1157 1159 WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1158 1160 WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1159 1161
1160 1162 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1161 1163
1162 1164 WijEijNij=numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1163 1165 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1164 1166
1165 1167 VxVy=numpy.array([[cA,cH],[cH,cB]])
1166 1168 VxVyResults=numpy.array([-cF,-cG])
1167 1169 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1168 1170
1169 1171 Vzon = Vy
1170 1172 Vmer = Vx
1171 1173 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1172 1174 Vang=numpy.arctan2(Vmer,Vzon)
1173 1175 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1174 1176 Vver=popt[1]
1175 1177 else:
1176 1178 Vver=numpy.NaN
1177 1179 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1178 1180
1179 1181
1180 1182 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1181 1183
1182 1184 class SpectralMoments(Operation):
1183 1185
1184 1186 '''
1185 1187 Function SpectralMoments()
1186 1188
1187 1189 Calculates moments (power, mean, standard deviation) and SNR of the signal
1188 1190
1189 1191 Type of dataIn: Spectra
1190 1192
1191 1193 Configuration Parameters:
1192 1194
1193 1195 dirCosx : Cosine director in X axis
1194 1196 dirCosy : Cosine director in Y axis
1195 1197
1196 1198 elevation :
1197 1199 azimuth :
1198 1200
1199 1201 Input:
1200 1202 channelList : simple channel list to select e.g. [2,3,7]
1201 1203 self.dataOut.data_pre : Spectral data
1202 1204 self.dataOut.abscissaList : List of frequencies
1203 1205 self.dataOut.noise : Noise level per channel
1204 1206
1205 1207 Affected:
1206 1208 self.dataOut.moments : Parameters per channel
1207 1209 self.dataOut.data_SNR : SNR per channel
1208 1210
1209 1211 '''
1210 1212
1211 1213 def run(self, dataOut):
1212 1214
1213 1215 #dataOut.data_pre = dataOut.data_pre[0]
1214 1216 data = dataOut.data_pre[0]
1215 1217 absc = dataOut.abscissaList[:-1]
1216 1218 noise = dataOut.noise
1217 1219 nChannel = data.shape[0]
1218 1220 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
1219 1221
1220 1222 for ind in range(nChannel):
1221 1223 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1222 1224
1223 1225 dataOut.moments = data_param[:,1:,:]
1224 1226 dataOut.data_SNR = data_param[:,0]
1225 1227 return
1226 1228
1227 1229 def __calculateMoments(self, oldspec, oldfreq, n0,
1228 1230 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
1229 1231
1230 1232 if (nicoh == None): nicoh = 1
1231 1233 if (graph == None): graph = 0
1232 1234 if (smooth == None): smooth = 0
1233 1235 elif (self.smooth < 3): smooth = 0
1234 1236
1235 1237 if (type1 == None): type1 = 0
1236 1238 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
1237 1239 if (snrth == None): snrth = -3
1238 1240 if (dc == None): dc = 0
1239 1241 if (aliasing == None): aliasing = 0
1240 1242 if (oldfd == None): oldfd = 0
1241 1243 if (wwauto == None): wwauto = 0
1242 1244
1243 1245 if (n0 < 1.e-20): n0 = 1.e-20
1244 1246
1245 1247 freq = oldfreq
1246 1248 vec_power = numpy.zeros(oldspec.shape[1])
1247 1249 vec_fd = numpy.zeros(oldspec.shape[1])
1248 1250 vec_w = numpy.zeros(oldspec.shape[1])
1249 1251 vec_snr = numpy.zeros(oldspec.shape[1])
1250 1252
1251 1253 oldspec = numpy.ma.masked_invalid(oldspec)
1252 1254
1253 1255 for ind in range(oldspec.shape[1]):
1254 1256
1255 1257 spec = oldspec[:,ind]
1256 1258 aux = spec*fwindow
1257 1259 max_spec = aux.max()
1258 1260 m = list(aux).index(max_spec)
1259 1261
1260 1262 #Smooth
1261 1263 if (smooth == 0): spec2 = spec
1262 1264 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1263 1265
1264 1266 # Calculo de Momentos
1265 1267 bb = spec2[range(m,spec2.size)]
1266 1268 bb = (bb<n0).nonzero()
1267 1269 bb = bb[0]
1268 1270
1269 1271 ss = spec2[range(0,m + 1)]
1270 1272 ss = (ss<n0).nonzero()
1271 1273 ss = ss[0]
1272 1274
1273 1275 if (bb.size == 0):
1274 1276 bb0 = spec.size - 1 - m
1275 1277 else:
1276 1278 bb0 = bb[0] - 1
1277 1279 if (bb0 < 0):
1278 1280 bb0 = 0
1279 1281
1280 1282 if (ss.size == 0): ss1 = 1
1281 1283 else: ss1 = max(ss) + 1
1282 1284
1283 1285 if (ss1 > m): ss1 = m
1284 1286
1285 1287 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
1286 1288 power = ( (spec2[valid] - n0) * fwindow[valid] ).sum()
1287 1289 fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power
1288 1290
1289 1291 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1290 1292 snr = (spec2.mean()-n0)/n0
1291 1293
1292 1294 if (snr < 1.e-20) :
1293 1295 snr = 1.e-20
1294 1296
1295 1297 vec_power[ind] = power
1296 1298 vec_fd[ind] = fd
1297 1299 vec_w[ind] = w
1298 1300 vec_snr[ind] = snr
1299 1301
1300 1302 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1301 1303 return moments
1302 1304
1303 1305 #------------------ Get SA Parameters --------------------------
1304 1306
1305 1307 def GetSAParameters(self):
1306 1308 #SA en frecuencia
1307 1309 pairslist = self.dataOut.groupList
1308 1310 num_pairs = len(pairslist)
1309 1311
1310 1312 vel = self.dataOut.abscissaList
1311 1313 spectra = self.dataOut.data_pre
1312 1314 cspectra = self.dataIn.data_cspc
1313 1315 delta_v = vel[1] - vel[0]
1314 1316
1315 1317 #Calculating the power spectrum
1316 1318 spc_pow = numpy.sum(spectra, 3)*delta_v
1317 1319 #Normalizing Spectra
1318 1320 norm_spectra = spectra/spc_pow
1319 1321 #Calculating the norm_spectra at peak
1320 1322 max_spectra = numpy.max(norm_spectra, 3)
1321 1323
1322 1324 #Normalizing Cross Spectra
1323 1325 norm_cspectra = numpy.zeros(cspectra.shape)
1324 1326
1325 1327 for i in range(num_chan):
1326 1328 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
1327 1329
1328 1330 max_cspectra = numpy.max(norm_cspectra,2)
1329 1331 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
1330 1332
1331 1333 for i in range(num_pairs):
1332 1334 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
1333 1335 #------------------- Get Lags ----------------------------------
1334 1336
1335 1337 class SALags(Operation):
1336 1338 '''
1337 1339 Function GetMoments()
1338 1340
1339 1341 Input:
1340 1342 self.dataOut.data_pre
1341 1343 self.dataOut.abscissaList
1342 1344 self.dataOut.noise
1343 1345 self.dataOut.normFactor
1344 1346 self.dataOut.data_SNR
1345 1347 self.dataOut.groupList
1346 1348 self.dataOut.nChannels
1347 1349
1348 1350 Affected:
1349 1351 self.dataOut.data_param
1350 1352
1351 1353 '''
1352 1354 def run(self, dataOut):
1353 1355 data_acf = dataOut.data_pre[0]
1354 1356 data_ccf = dataOut.data_pre[1]
1355 1357 normFactor_acf = dataOut.normFactor[0]
1356 1358 normFactor_ccf = dataOut.normFactor[1]
1357 1359 pairs_acf = dataOut.groupList[0]
1358 1360 pairs_ccf = dataOut.groupList[1]
1359 1361
1360 1362 nHeights = dataOut.nHeights
1361 1363 absc = dataOut.abscissaList
1362 1364 noise = dataOut.noise
1363 1365 SNR = dataOut.data_SNR
1364 1366 nChannels = dataOut.nChannels
1365 1367 # pairsList = dataOut.groupList
1366 1368 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1367 1369
1368 1370 for l in range(len(pairs_acf)):
1369 1371 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
1370 1372
1371 1373 for l in range(len(pairs_ccf)):
1372 1374 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
1373 1375
1374 1376 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
1375 1377 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
1376 1378 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
1377 1379 return
1378 1380
1379 1381 # def __getPairsAutoCorr(self, pairsList, nChannels):
1380 1382 #
1381 1383 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1382 1384 #
1383 1385 # for l in range(len(pairsList)):
1384 1386 # firstChannel = pairsList[l][0]
1385 1387 # secondChannel = pairsList[l][1]
1386 1388 #
1387 1389 # #Obteniendo pares de Autocorrelacion
1388 1390 # if firstChannel == secondChannel:
1389 1391 # pairsAutoCorr[firstChannel] = int(l)
1390 1392 #
1391 1393 # pairsAutoCorr = pairsAutoCorr.astype(int)
1392 1394 #
1393 1395 # pairsCrossCorr = range(len(pairsList))
1394 1396 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1395 1397 #
1396 1398 # return pairsAutoCorr, pairsCrossCorr
1397 1399
1398 1400 def __calculateTaus(self, data_acf, data_ccf, lagRange):
1399 1401
1400 1402 lag0 = data_acf.shape[1]/2
1401 1403 #Funcion de Autocorrelacion
1402 1404 mean_acf = stats.nanmean(data_acf, axis = 0)
1403 1405
1404 1406 #Obtencion Indice de TauCross
1405 1407 ind_ccf = data_ccf.argmax(axis = 1)
1406 1408 #Obtencion Indice de TauAuto
1407 1409 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
1408 1410 ccf_lag0 = data_ccf[:,lag0,:]
1409 1411
1410 1412 for i in range(ccf_lag0.shape[0]):
1411 1413 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
1412 1414
1413 1415 #Obtencion de TauCross y TauAuto
1414 1416 tau_ccf = lagRange[ind_ccf]
1415 1417 tau_acf = lagRange[ind_acf]
1416 1418
1417 1419 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
1418 1420
1419 1421 tau_ccf[Nan1,Nan2] = numpy.nan
1420 1422 tau_acf[Nan1,Nan2] = numpy.nan
1421 1423 tau = numpy.vstack((tau_ccf,tau_acf))
1422 1424
1423 1425 return tau
1424 1426
1425 1427 def __calculateLag1Phase(self, data, lagTRange):
1426 1428 data1 = stats.nanmean(data, axis = 0)
1427 1429 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
1428 1430
1429 1431 phase = numpy.angle(data1[lag1,:])
1430 1432
1431 1433 return phase
1432 1434
1433 1435 class SpectralFitting(Operation):
1434 1436 '''
1435 1437 Function GetMoments()
1436 1438
1437 1439 Input:
1438 1440 Output:
1439 1441 Variables modified:
1440 1442 '''
1441 1443
1442 1444 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1443 1445
1444 1446
1445 1447 if path != None:
1446 1448 sys.path.append(path)
1447 1449 self.dataOut.library = importlib.import_module(file)
1448 1450
1449 1451 #To be inserted as a parameter
1450 1452 groupArray = numpy.array(groupList)
1451 1453 # groupArray = numpy.array([[0,1],[2,3]])
1452 1454 self.dataOut.groupList = groupArray
1453 1455
1454 1456 nGroups = groupArray.shape[0]
1455 1457 nChannels = self.dataIn.nChannels
1456 1458 nHeights=self.dataIn.heightList.size
1457 1459
1458 1460 #Parameters Array
1459 1461 self.dataOut.data_param = None
1460 1462
1461 1463 #Set constants
1462 1464 constants = self.dataOut.library.setConstants(self.dataIn)
1463 1465 self.dataOut.constants = constants
1464 1466 M = self.dataIn.normFactor
1465 1467 N = self.dataIn.nFFTPoints
1466 1468 ippSeconds = self.dataIn.ippSeconds
1467 1469 K = self.dataIn.nIncohInt
1468 1470 pairsArray = numpy.array(self.dataIn.pairsList)
1469 1471
1470 1472 #List of possible combinations
1471 1473 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1472 1474 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1473 1475
1474 1476 if getSNR:
1475 1477 listChannels = groupArray.reshape((groupArray.size))
1476 1478 listChannels.sort()
1477 1479 noise = self.dataIn.getNoise()
1478 1480 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1479 1481
1480 1482 for i in range(nGroups):
1481 1483 coord = groupArray[i,:]
1482 1484
1483 1485 #Input data array
1484 1486 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1485 1487 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1486 1488
1487 1489 #Cross Spectra data array for Covariance Matrixes
1488 1490 ind = 0
1489 1491 for pairs in listComb:
1490 1492 pairsSel = numpy.array([coord[x],coord[y]])
1491 1493 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1492 1494 ind += 1
1493 1495 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1494 1496 dataCross = dataCross**2/K
1495 1497
1496 1498 for h in range(nHeights):
1497 1499
1498 1500 #Input
1499 1501 d = data[:,h]
1500 1502
1501 1503 #Covariance Matrix
1502 1504 D = numpy.diag(d**2/K)
1503 1505 ind = 0
1504 1506 for pairs in listComb:
1505 1507 #Coordinates in Covariance Matrix
1506 1508 x = pairs[0]
1507 1509 y = pairs[1]
1508 1510 #Channel Index
1509 1511 S12 = dataCross[ind,:,h]
1510 1512 D12 = numpy.diag(S12)
1511 1513 #Completing Covariance Matrix with Cross Spectras
1512 1514 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1513 1515 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1514 1516 ind += 1
1515 1517 Dinv=numpy.linalg.inv(D)
1516 1518 L=numpy.linalg.cholesky(Dinv)
1517 1519 LT=L.T
1518 1520
1519 1521 dp = numpy.dot(LT,d)
1520 1522
1521 1523 #Initial values
1522 1524 data_spc = self.dataIn.data_spc[coord,:,h]
1523 1525
1524 1526 if (h>0)and(error1[3]<5):
1525 1527 p0 = self.dataOut.data_param[i,:,h-1]
1526 1528 else:
1527 1529 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1528 1530
1529 1531 try:
1530 1532 #Least Squares
1531 1533 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1532 1534 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1533 1535 #Chi square error
1534 1536 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1535 1537 #Error with Jacobian
1536 1538 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1537 1539 except:
1538 1540 minp = p0*numpy.nan
1539 1541 error0 = numpy.nan
1540 1542 error1 = p0*numpy.nan
1541 1543
1542 1544 #Save
1543 1545 if self.dataOut.data_param == None:
1544 1546 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1545 1547 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1546 1548
1547 1549 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1548 1550 self.dataOut.data_param[i,:,h] = minp
1549 1551 return
1550 1552
1551 1553 def __residFunction(self, p, dp, LT, constants):
1552 1554
1553 1555 fm = self.dataOut.library.modelFunction(p, constants)
1554 1556 fmp=numpy.dot(LT,fm)
1555 1557
1556 1558 return dp-fmp
1557 1559
1558 1560 def __getSNR(self, z, noise):
1559 1561
1560 1562 avg = numpy.average(z, axis=1)
1561 1563 SNR = (avg.T-noise)/noise
1562 1564 SNR = SNR.T
1563 1565 return SNR
1564 1566
1565 1567 def __chisq(p,chindex,hindex):
1566 1568 #similar to Resid but calculates CHI**2
1567 1569 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1568 1570 dp=numpy.dot(LT,d)
1569 1571 fmp=numpy.dot(LT,fm)
1570 1572 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1571 1573 return chisq
1572 1574
1573 1575 class WindProfiler(Operation):
1574 1576
1575 1577 __isConfig = False
1576 1578
1577 1579 __initime = None
1578 1580 __lastdatatime = None
1579 1581 __integrationtime = None
1580 1582
1581 1583 __buffer = None
1582 1584
1583 1585 __dataReady = False
1584 1586
1585 1587 __firstdata = None
1586 1588
1587 1589 n = None
1588 1590
1589 1591 def __init__(self, **kwargs):
1590 1592 Operation.__init__(self, **kwargs)
1591 1593
1592 1594 def __calculateCosDir(self, elev, azim):
1593 1595 zen = (90 - elev)*numpy.pi/180
1594 1596 azim = azim*numpy.pi/180
1595 1597 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1596 1598 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1597 1599
1598 1600 signX = numpy.sign(numpy.cos(azim))
1599 1601 signY = numpy.sign(numpy.sin(azim))
1600 1602
1601 1603 cosDirX = numpy.copysign(cosDirX, signX)
1602 1604 cosDirY = numpy.copysign(cosDirY, signY)
1603 1605 return cosDirX, cosDirY
1604 1606
1605 1607 def __calculateAngles(self, theta_x, theta_y, azimuth):
1606 1608
1607 1609 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1608 1610 zenith_arr = numpy.arccos(dir_cosw)
1609 1611 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1610 1612
1611 1613 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1612 1614 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1613 1615
1614 1616 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1615 1617
1616 1618 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1617 1619
1618 1620 #
1619 1621 if horOnly:
1620 1622 A = numpy.c_[dir_cosu,dir_cosv]
1621 1623 else:
1622 1624 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1623 1625 A = numpy.asmatrix(A)
1624 1626 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1625 1627
1626 1628 return A1
1627 1629
1628 1630 def __correctValues(self, heiRang, phi, velRadial, SNR):
1629 1631 listPhi = phi.tolist()
1630 1632 maxid = listPhi.index(max(listPhi))
1631 1633 minid = listPhi.index(min(listPhi))
1632 1634
1633 1635 rango = range(len(phi))
1634 1636 # rango = numpy.delete(rango,maxid)
1635 1637
1636 1638 heiRang1 = heiRang*math.cos(phi[maxid])
1637 1639 heiRangAux = heiRang*math.cos(phi[minid])
1638 1640 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1639 1641 heiRang1 = numpy.delete(heiRang1,indOut)
1640 1642
1641 1643 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1642 1644 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1643 1645
1644 1646 for i in rango:
1645 1647 x = heiRang*math.cos(phi[i])
1646 1648 y1 = velRadial[i,:]
1647 1649 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1648 1650
1649 1651 x1 = heiRang1
1650 1652 y11 = f1(x1)
1651 1653
1652 1654 y2 = SNR[i,:]
1653 1655 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1654 1656 y21 = f2(x1)
1655 1657
1656 1658 velRadial1[i,:] = y11
1657 1659 SNR1[i,:] = y21
1658 1660
1659 1661 return heiRang1, velRadial1, SNR1
1660 1662
1661 1663 def __calculateVelUVW(self, A, velRadial):
1662 1664
1663 1665 #Operacion Matricial
1664 1666 # velUVW = numpy.zeros((velRadial.shape[1],3))
1665 1667 # for ind in range(velRadial.shape[1]):
1666 1668 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1667 1669 # velUVW = velUVW.transpose()
1668 1670 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1669 1671 velUVW[:,:] = numpy.dot(A,velRadial)
1670 1672
1671 1673
1672 1674 return velUVW
1673 1675
1674 1676 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1675 1677
1676 1678 def techniqueDBS(self, kwargs):
1677 1679 """
1678 1680 Function that implements Doppler Beam Swinging (DBS) technique.
1679 1681
1680 1682 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1681 1683 Direction correction (if necessary), Ranges and SNR
1682 1684
1683 1685 Output: Winds estimation (Zonal, Meridional and Vertical)
1684 1686
1685 1687 Parameters affected: Winds, height range, SNR
1686 1688 """
1687 1689 velRadial0 = kwargs['velRadial']
1688 1690 heiRang = kwargs['heightList']
1689 1691 SNR0 = kwargs['SNR']
1690 1692
1691 1693 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1692 1694 theta_x = numpy.array(kwargs['dirCosx'])
1693 1695 theta_y = numpy.array(kwargs['dirCosy'])
1694 1696 else:
1695 1697 elev = numpy.array(kwargs['elevation'])
1696 1698 azim = numpy.array(kwargs['azimuth'])
1697 1699 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1698 1700 azimuth = kwargs['correctAzimuth']
1699 1701 if kwargs.has_key('horizontalOnly'):
1700 1702 horizontalOnly = kwargs['horizontalOnly']
1701 1703 else: horizontalOnly = False
1702 1704 if kwargs.has_key('correctFactor'):
1703 1705 correctFactor = kwargs['correctFactor']
1704 1706 else: correctFactor = 1
1705 1707 if kwargs.has_key('channelList'):
1706 1708 channelList = kwargs['channelList']
1707 1709 if len(channelList) == 2:
1708 1710 horizontalOnly = True
1709 1711 arrayChannel = numpy.array(channelList)
1710 1712 param = param[arrayChannel,:,:]
1711 1713 theta_x = theta_x[arrayChannel]
1712 1714 theta_y = theta_y[arrayChannel]
1713 1715
1714 1716 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1715 1717 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
1716 1718 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1717 1719
1718 1720 #Calculo de Componentes de la velocidad con DBS
1719 1721 winds = self.__calculateVelUVW(A,velRadial1)
1720 1722
1721 1723 return winds, heiRang1, SNR1
1722 1724
1723 1725 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
1724 1726
1725 1727 nPairs = len(pairs_ccf)
1726 1728 posx = numpy.asarray(posx)
1727 1729 posy = numpy.asarray(posy)
1728 1730
1729 1731 #Rotacion Inversa para alinear con el azimuth
1730 1732 if azimuth!= None:
1731 1733 azimuth = azimuth*math.pi/180
1732 1734 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1733 1735 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1734 1736 else:
1735 1737 posx1 = posx
1736 1738 posy1 = posy
1737 1739
1738 1740 #Calculo de Distancias
1739 1741 distx = numpy.zeros(nPairs)
1740 1742 disty = numpy.zeros(nPairs)
1741 1743 dist = numpy.zeros(nPairs)
1742 1744 ang = numpy.zeros(nPairs)
1743 1745
1744 1746 for i in range(nPairs):
1745 1747 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
1746 1748 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
1747 1749 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1748 1750 ang[i] = numpy.arctan2(disty[i],distx[i])
1749 1751
1750 1752 return distx, disty, dist, ang
1751 1753 #Calculo de Matrices
1752 1754 # nPairs = len(pairs)
1753 1755 # ang1 = numpy.zeros((nPairs, 2, 1))
1754 1756 # dist1 = numpy.zeros((nPairs, 2, 1))
1755 1757 #
1756 1758 # for j in range(nPairs):
1757 1759 # dist1[j,0,0] = dist[pairs[j][0]]
1758 1760 # dist1[j,1,0] = dist[pairs[j][1]]
1759 1761 # ang1[j,0,0] = ang[pairs[j][0]]
1760 1762 # ang1[j,1,0] = ang[pairs[j][1]]
1761 1763 #
1762 1764 # return distx,disty, dist1,ang1
1763 1765
1764 1766
1765 1767 def __calculateVelVer(self, phase, lagTRange, _lambda):
1766 1768
1767 1769 Ts = lagTRange[1] - lagTRange[0]
1768 1770 velW = -_lambda*phase/(4*math.pi*Ts)
1769 1771
1770 1772 return velW
1771 1773
1772 1774 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1773 1775 nPairs = tau1.shape[0]
1774 1776 nHeights = tau1.shape[1]
1775 1777 vel = numpy.zeros((nPairs,3,nHeights))
1776 1778 dist1 = numpy.reshape(dist, (dist.size,1))
1777 1779
1778 1780 angCos = numpy.cos(ang)
1779 1781 angSin = numpy.sin(ang)
1780 1782
1781 1783 vel0 = dist1*tau1/(2*tau2**2)
1782 1784 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1783 1785 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1784 1786
1785 1787 ind = numpy.where(numpy.isinf(vel))
1786 1788 vel[ind] = numpy.nan
1787 1789
1788 1790 return vel
1789 1791
1790 1792 # def __getPairsAutoCorr(self, pairsList, nChannels):
1791 1793 #
1792 1794 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1793 1795 #
1794 1796 # for l in range(len(pairsList)):
1795 1797 # firstChannel = pairsList[l][0]
1796 1798 # secondChannel = pairsList[l][1]
1797 1799 #
1798 1800 # #Obteniendo pares de Autocorrelacion
1799 1801 # if firstChannel == secondChannel:
1800 1802 # pairsAutoCorr[firstChannel] = int(l)
1801 1803 #
1802 1804 # pairsAutoCorr = pairsAutoCorr.astype(int)
1803 1805 #
1804 1806 # pairsCrossCorr = range(len(pairsList))
1805 1807 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1806 1808 #
1807 1809 # return pairsAutoCorr, pairsCrossCorr
1808 1810
1809 1811 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1810 1812 def techniqueSA(self, kwargs):
1811 1813
1812 1814 """
1813 1815 Function that implements Spaced Antenna (SA) technique.
1814 1816
1815 1817 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1816 1818 Direction correction (if necessary), Ranges and SNR
1817 1819
1818 1820 Output: Winds estimation (Zonal, Meridional and Vertical)
1819 1821
1820 1822 Parameters affected: Winds
1821 1823 """
1822 1824 position_x = kwargs['positionX']
1823 1825 position_y = kwargs['positionY']
1824 1826 azimuth = kwargs['azimuth']
1825 1827
1826 1828 if kwargs.has_key('correctFactor'):
1827 1829 correctFactor = kwargs['correctFactor']
1828 1830 else:
1829 1831 correctFactor = 1
1830 1832
1831 1833 groupList = kwargs['groupList']
1832 1834 pairs_ccf = groupList[1]
1833 1835 tau = kwargs['tau']
1834 1836 _lambda = kwargs['_lambda']
1835 1837
1836 1838 #Cross Correlation pairs obtained
1837 1839 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
1838 1840 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1839 1841 # pairsSelArray = numpy.array(pairsSelected)
1840 1842 # pairs = []
1841 1843 #
1842 1844 # #Wind estimation pairs obtained
1843 1845 # for i in range(pairsSelArray.shape[0]/2):
1844 1846 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1845 1847 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1846 1848 # pairs.append((ind1,ind2))
1847 1849
1848 1850 indtau = tau.shape[0]/2
1849 1851 tau1 = tau[:indtau,:]
1850 1852 tau2 = tau[indtau:-1,:]
1851 1853 # tau1 = tau1[pairs,:]
1852 1854 # tau2 = tau2[pairs,:]
1853 1855 phase1 = tau[-1,:]
1854 1856
1855 1857 #---------------------------------------------------------------------
1856 1858 #Metodo Directo
1857 1859 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
1858 1860 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1859 1861 winds = stats.nanmean(winds, axis=0)
1860 1862 #---------------------------------------------------------------------
1861 1863 #Metodo General
1862 1864 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1863 1865 # #Calculo Coeficientes de Funcion de Correlacion
1864 1866 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1865 1867 # #Calculo de Velocidades
1866 1868 # winds = self.calculateVelUV(F,G,A,B,H)
1867 1869
1868 1870 #---------------------------------------------------------------------
1869 1871 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1870 1872 winds = correctFactor*winds
1871 1873 return winds
1872 1874
1873 1875 def __checkTime(self, currentTime, paramInterval, outputInterval):
1874 1876
1875 1877 dataTime = currentTime + paramInterval
1876 1878 deltaTime = dataTime - self.__initime
1877 1879
1878 1880 if deltaTime >= outputInterval or deltaTime < 0:
1879 1881 self.__dataReady = True
1880 1882 return
1881 1883
1882 1884 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1883 1885 '''
1884 1886 Function that implements winds estimation technique with detected meteors.
1885 1887
1886 1888 Input: Detected meteors, Minimum meteor quantity to wind estimation
1887 1889
1888 1890 Output: Winds estimation (Zonal and Meridional)
1889 1891
1890 1892 Parameters affected: Winds
1891 1893 '''
1892 1894 # print arrayMeteor.shape
1893 1895 #Settings
1894 1896 nInt = (heightMax - heightMin)/2
1895 1897 # print nInt
1896 1898 nInt = int(nInt)
1897 1899 # print nInt
1898 1900 winds = numpy.zeros((2,nInt))*numpy.nan
1899 1901
1900 1902 #Filter errors
1901 1903 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
1902 1904 finalMeteor = arrayMeteor[error,:]
1903 1905
1904 1906 #Meteor Histogram
1905 1907 finalHeights = finalMeteor[:,2]
1906 1908 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1907 1909 nMeteorsPerI = hist[0]
1908 1910 heightPerI = hist[1]
1909 1911
1910 1912 #Sort of meteors
1911 1913 indSort = finalHeights.argsort()
1912 1914 finalMeteor2 = finalMeteor[indSort,:]
1913 1915
1914 1916 # Calculating winds
1915 1917 ind1 = 0
1916 1918 ind2 = 0
1917 1919
1918 1920 for i in range(nInt):
1919 1921 nMet = nMeteorsPerI[i]
1920 1922 ind1 = ind2
1921 1923 ind2 = ind1 + nMet
1922 1924
1923 1925 meteorAux = finalMeteor2[ind1:ind2,:]
1924 1926
1925 1927 if meteorAux.shape[0] >= meteorThresh:
1926 1928 vel = meteorAux[:, 6]
1927 1929 zen = meteorAux[:, 4]*numpy.pi/180
1928 1930 azim = meteorAux[:, 3]*numpy.pi/180
1929 1931
1930 1932 n = numpy.cos(zen)
1931 1933 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1932 1934 # l = m*numpy.tan(azim)
1933 1935 l = numpy.sin(zen)*numpy.sin(azim)
1934 1936 m = numpy.sin(zen)*numpy.cos(azim)
1935 1937
1936 1938 A = numpy.vstack((l, m)).transpose()
1937 1939 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1938 1940 windsAux = numpy.dot(A1, vel)
1939 1941
1940 1942 winds[0,i] = windsAux[0]
1941 1943 winds[1,i] = windsAux[1]
1942 1944
1943 1945 return winds, heightPerI[:-1]
1944 1946
1945 1947 def techniqueNSM_SA(self, **kwargs):
1946 1948 metArray = kwargs['metArray']
1947 1949 heightList = kwargs['heightList']
1948 1950 timeList = kwargs['timeList']
1949 1951
1950 1952 rx_location = kwargs['rx_location']
1951 1953 groupList = kwargs['groupList']
1952 1954 azimuth = kwargs['azimuth']
1953 1955 dfactor = kwargs['dfactor']
1954 1956 k = kwargs['k']
1955 1957
1956 1958 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
1957 1959 d = dist*dfactor
1958 1960 #Phase calculation
1959 1961 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
1960 1962
1961 1963 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
1962 1964
1963 1965 velEst = numpy.zeros((heightList.size,2))*numpy.nan
1964 1966 azimuth1 = azimuth1*numpy.pi/180
1965 1967
1966 1968 for i in range(heightList.size):
1967 1969 h = heightList[i]
1968 1970 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
1969 1971 metHeight = metArray1[indH,:]
1970 1972 if metHeight.shape[0] >= 2:
1971 1973 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
1972 1974 iazim = metHeight[:,1].astype(int)
1973 1975 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
1974 1976 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
1975 1977 A = numpy.asmatrix(A)
1976 1978 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
1977 1979 velHor = numpy.dot(A1,velAux)
1978 1980
1979 1981 velEst[i,:] = numpy.squeeze(velHor)
1980 1982 return velEst
1981 1983
1982 1984 def __getPhaseSlope(self, metArray, heightList, timeList):
1983 1985 meteorList = []
1984 1986 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
1985 1987 #Putting back together the meteor matrix
1986 1988 utctime = metArray[:,0]
1987 1989 uniqueTime = numpy.unique(utctime)
1988 1990
1989 1991 phaseDerThresh = 0.5
1990 1992 ippSeconds = timeList[1] - timeList[0]
1991 1993 sec = numpy.where(timeList>1)[0][0]
1992 1994 nPairs = metArray.shape[1] - 6
1993 1995 nHeights = len(heightList)
1994 1996
1995 1997 for t in uniqueTime:
1996 1998 metArray1 = metArray[utctime==t,:]
1997 1999 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
1998 2000 tmet = metArray1[:,1].astype(int)
1999 2001 hmet = metArray1[:,2].astype(int)
2000 2002
2001 2003 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
2002 2004 metPhase[:,:] = numpy.nan
2003 2005 metPhase[:,hmet,tmet] = metArray1[:,6:].T
2004 2006
2005 2007 #Delete short trails
2006 2008 metBool = ~numpy.isnan(metPhase[0,:,:])
2007 2009 heightVect = numpy.sum(metBool, axis = 1)
2008 2010 metBool[heightVect<sec,:] = False
2009 2011 metPhase[:,heightVect<sec,:] = numpy.nan
2010 2012
2011 2013 #Derivative
2012 2014 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
2013 2015 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
2014 2016 metPhase[phDerAux] = numpy.nan
2015 2017
2016 2018 #--------------------------METEOR DETECTION -----------------------------------------
2017 2019 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
2018 2020
2019 2021 for p in numpy.arange(nPairs):
2020 2022 phase = metPhase[p,:,:]
2021 2023 phDer = metDer[p,:,:]
2022 2024
2023 2025 for h in indMet:
2024 2026 height = heightList[h]
2025 2027 phase1 = phase[h,:] #82
2026 2028 phDer1 = phDer[h,:]
2027 2029
2028 2030 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
2029 2031
2030 2032 indValid = numpy.where(~numpy.isnan(phase1))[0]
2031 2033 initMet = indValid[0]
2032 2034 endMet = 0
2033 2035
2034 2036 for i in range(len(indValid)-1):
2035 2037
2036 2038 #Time difference
2037 2039 inow = indValid[i]
2038 2040 inext = indValid[i+1]
2039 2041 idiff = inext - inow
2040 2042 #Phase difference
2041 2043 phDiff = numpy.abs(phase1[inext] - phase1[inow])
2042 2044
2043 2045 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
2044 2046 sizeTrail = inow - initMet + 1
2045 2047 if sizeTrail>3*sec: #Too short meteors
2046 2048 x = numpy.arange(initMet,inow+1)*ippSeconds
2047 2049 y = phase1[initMet:inow+1]
2048 2050 ynnan = ~numpy.isnan(y)
2049 2051 x = x[ynnan]
2050 2052 y = y[ynnan]
2051 2053 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
2052 2054 ylin = x*slope + intercept
2053 2055 rsq = r_value**2
2054 2056 if rsq > 0.5:
2055 2057 vel = slope#*height*1000/(k*d)
2056 2058 estAux = numpy.array([utctime,p,height, vel, rsq])
2057 2059 meteorList.append(estAux)
2058 2060 initMet = inext
2059 2061 metArray2 = numpy.array(meteorList)
2060 2062
2061 2063 return metArray2
2062 2064
2063 2065 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
2064 2066
2065 2067 azimuth1 = numpy.zeros(len(pairslist))
2066 2068 dist = numpy.zeros(len(pairslist))
2067 2069
2068 2070 for i in range(len(rx_location)):
2069 2071 ch0 = pairslist[i][0]
2070 2072 ch1 = pairslist[i][1]
2071 2073
2072 2074 diffX = rx_location[ch0][0] - rx_location[ch1][0]
2073 2075 diffY = rx_location[ch0][1] - rx_location[ch1][1]
2074 2076 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
2075 2077 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
2076 2078
2077 2079 azimuth1 -= azimuth0
2078 2080 return azimuth1, dist
2079 2081
2080 2082 def techniqueNSM_DBS(self, **kwargs):
2081 2083 metArray = kwargs['metArray']
2082 2084 heightList = kwargs['heightList']
2083 2085 timeList = kwargs['timeList']
2084 2086 zenithList = kwargs['zenithList']
2085 2087 nChan = numpy.max(cmet) + 1
2086 2088 nHeights = len(heightList)
2087 2089
2088 2090 utctime = metArray[:,0]
2089 2091 cmet = metArray[:,1]
2090 2092 hmet = metArray1[:,3].astype(int)
2091 2093 h1met = heightList[hmet]*zenithList[cmet]
2092 2094 vmet = metArray1[:,5]
2093 2095
2094 2096 for i in range(nHeights - 1):
2095 2097 hmin = heightList[i]
2096 2098 hmax = heightList[i + 1]
2097 2099
2098 2100 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
2099 2101
2100 2102
2101 2103
2102 2104 return data_output
2103 2105
2104 2106 def run(self, dataOut, technique, positionY, positionX, azimuth, **kwargs):
2105 2107
2106 2108 param = dataOut.data_param
2107 2109 if dataOut.abscissaList != None:
2108 2110 absc = dataOut.abscissaList[:-1]
2109 2111 noise = dataOut.noise
2110 2112 heightList = dataOut.heightList
2111 2113 SNR = dataOut.data_SNR
2112 2114
2113 2115 if technique == 'DBS':
2114 2116
2115 2117 kwargs['velRadial'] = param[:,1,:] #Radial velocity
2116 2118 kwargs['heightList'] = heightList
2117 2119 kwargs['SNR'] = SNR
2118 2120
2119 2121 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
2120 2122 dataOut.utctimeInit = dataOut.utctime
2121 2123 dataOut.outputInterval = dataOut.paramInterval
2122 2124
2123 2125 elif technique == 'SA':
2124 2126
2125 2127 #Parameters
2126 2128 # position_x = kwargs['positionX']
2127 2129 # position_y = kwargs['positionY']
2128 2130 # azimuth = kwargs['azimuth']
2129 2131 #
2130 2132 # if kwargs.has_key('crosspairsList'):
2131 2133 # pairs = kwargs['crosspairsList']
2132 2134 # else:
2133 2135 # pairs = None
2134 2136 #
2135 2137 # if kwargs.has_key('correctFactor'):
2136 2138 # correctFactor = kwargs['correctFactor']
2137 2139 # else:
2138 2140 # correctFactor = 1
2139 2141
2140 2142 # tau = dataOut.data_param
2141 2143 # _lambda = dataOut.C/dataOut.frequency
2142 2144 # pairsList = dataOut.groupList
2143 2145 # nChannels = dataOut.nChannels
2144 2146
2145 2147 kwargs['groupList'] = dataOut.groupList
2146 2148 kwargs['tau'] = dataOut.data_param
2147 2149 kwargs['_lambda'] = dataOut.C/dataOut.frequency
2148 2150 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
2149 2151 dataOut.data_output = self.techniqueSA(kwargs)
2150 2152 dataOut.utctimeInit = dataOut.utctime
2151 2153 dataOut.outputInterval = dataOut.timeInterval
2152 2154
2153 2155 elif technique == 'Meteors':
2154 2156 dataOut.flagNoData = True
2155 2157 self.__dataReady = False
2156 2158
2157 2159 if kwargs.has_key('nHours'):
2158 2160 nHours = kwargs['nHours']
2159 2161 else:
2160 2162 nHours = 1
2161 2163
2162 2164 if kwargs.has_key('meteorsPerBin'):
2163 2165 meteorThresh = kwargs['meteorsPerBin']
2164 2166 else:
2165 2167 meteorThresh = 6
2166 2168
2167 2169 if kwargs.has_key('hmin'):
2168 2170 hmin = kwargs['hmin']
2169 2171 else: hmin = 70
2170 2172 if kwargs.has_key('hmax'):
2171 2173 hmax = kwargs['hmax']
2172 2174 else: hmax = 110
2173 2175
2174 2176 dataOut.outputInterval = nHours*3600
2175 2177
2176 2178 if self.__isConfig == False:
2177 2179 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2178 2180 #Get Initial LTC time
2179 2181 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2180 2182 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2181 2183
2182 2184 self.__isConfig = True
2183 2185
2184 2186 if self.__buffer == None:
2185 2187 self.__buffer = dataOut.data_param
2186 2188 self.__firstdata = copy.copy(dataOut)
2187 2189
2188 2190 else:
2189 2191 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2190 2192
2191 2193 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2192 2194
2193 2195 if self.__dataReady:
2194 2196 dataOut.utctimeInit = self.__initime
2195 2197
2196 2198 self.__initime += dataOut.outputInterval #to erase time offset
2197 2199
2198 2200 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
2199 2201 dataOut.flagNoData = False
2200 2202 self.__buffer = None
2201 2203
2202 2204 elif technique == 'Meteors1':
2203 2205 dataOut.flagNoData = True
2204 2206 self.__dataReady = False
2205 2207
2206 2208 if kwargs.has_key('nMins'):
2207 2209 nMins = kwargs['nMins']
2208 2210 else: nMins = 20
2209 2211 if kwargs.has_key('rx_location'):
2210 2212 rx_location = kwargs['rx_location']
2211 2213 else: rx_location = [(0,1),(1,1),(1,0)]
2212 2214 if kwargs.has_key('azimuth'):
2213 2215 azimuth = kwargs['azimuth']
2214 2216 else: azimuth = 51
2215 2217 if kwargs.has_key('dfactor'):
2216 2218 dfactor = kwargs['dfactor']
2217 2219 if kwargs.has_key('mode'):
2218 2220 mode = kwargs['mode']
2219 2221 else: mode = 'SA'
2220 2222
2221 2223 #Borrar luego esto
2222 2224 if dataOut.groupList == None:
2223 2225 dataOut.groupList = [(0,1),(0,2),(1,2)]
2224 2226 groupList = dataOut.groupList
2225 2227 C = 3e8
2226 2228 freq = 50e6
2227 2229 lamb = C/freq
2228 2230 k = 2*numpy.pi/lamb
2229 2231
2230 2232 timeList = dataOut.abscissaList
2231 2233 heightList = dataOut.heightList
2232 2234
2233 2235 if self.__isConfig == False:
2234 2236 dataOut.outputInterval = nMins*60
2235 2237 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2236 2238 #Get Initial LTC time
2237 2239 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2238 2240 minuteAux = initime.minute
2239 2241 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
2240 2242 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2241 2243
2242 2244 self.__isConfig = True
2243 2245
2244 2246 if self.__buffer == None:
2245 2247 self.__buffer = dataOut.data_param
2246 2248 self.__firstdata = copy.copy(dataOut)
2247 2249
2248 2250 else:
2249 2251 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2250 2252
2251 2253 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2252 2254
2253 2255 if self.__dataReady:
2254 2256 dataOut.utctimeInit = self.__initime
2255 2257 self.__initime += dataOut.outputInterval #to erase time offset
2256 2258
2257 2259 metArray = self.__buffer
2258 2260 if mode == 'SA':
2259 2261 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
2260 2262 elif mode == 'DBS':
2261 2263 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
2262 2264 dataOut.data_output = dataOut.data_output.T
2263 2265 dataOut.flagNoData = False
2264 2266 self.__buffer = None
2265 2267
2266 2268 return
2267 2269
2268 2270 class EWDriftsEstimation(Operation):
2269 2271
2270 2272 def __init__(self):
2271 2273 Operation.__init__(self)
2272 2274
2273 2275 def __correctValues(self, heiRang, phi, velRadial, SNR):
2274 2276 listPhi = phi.tolist()
2275 2277 maxid = listPhi.index(max(listPhi))
2276 2278 minid = listPhi.index(min(listPhi))
2277 2279
2278 2280 rango = range(len(phi))
2279 2281 # rango = numpy.delete(rango,maxid)
2280 2282
2281 2283 heiRang1 = heiRang*math.cos(phi[maxid])
2282 2284 heiRangAux = heiRang*math.cos(phi[minid])
2283 2285 indOut = (heiRang1 < heiRangAux[0]).nonzero()
2284 2286 heiRang1 = numpy.delete(heiRang1,indOut)
2285 2287
2286 2288 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2287 2289 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2288 2290
2289 2291 for i in rango:
2290 2292 x = heiRang*math.cos(phi[i])
2291 2293 y1 = velRadial[i,:]
2292 2294 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2293 2295
2294 2296 x1 = heiRang1
2295 2297 y11 = f1(x1)
2296 2298
2297 2299 y2 = SNR[i,:]
2298 2300 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2299 2301 y21 = f2(x1)
2300 2302
2301 2303 velRadial1[i,:] = y11
2302 2304 SNR1[i,:] = y21
2303 2305
2304 2306 return heiRang1, velRadial1, SNR1
2305 2307
2306 2308 def run(self, dataOut, zenith, zenithCorrection):
2307 2309 heiRang = dataOut.heightList
2308 2310 velRadial = dataOut.data_param[:,3,:]
2309 2311 SNR = dataOut.data_SNR
2310 2312
2311 2313 zenith = numpy.array(zenith)
2312 2314 zenith -= zenithCorrection
2313 2315 zenith *= numpy.pi/180
2314 2316
2315 2317 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2316 2318
2317 2319 alp = zenith[0]
2318 2320 bet = zenith[1]
2319 2321
2320 2322 w_w = velRadial1[0,:]
2321 2323 w_e = velRadial1[1,:]
2322 2324
2323 2325 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
2324 2326 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
2325 2327
2326 2328 winds = numpy.vstack((u,w))
2327 2329
2328 2330 dataOut.heightList = heiRang1
2329 2331 dataOut.data_output = winds
2330 2332 dataOut.data_SNR = SNR1
2331 2333
2332 2334 dataOut.utctimeInit = dataOut.utctime
2333 2335 dataOut.outputInterval = dataOut.timeInterval
2334 2336 return
2335 2337
2336 2338 #--------------- Non Specular Meteor ----------------
2337 2339
2338 2340 class NonSpecularMeteorDetection(Operation):
2339 2341
2340 2342 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
2341 2343 data_acf = self.dataOut.data_pre[0]
2342 2344 data_ccf = self.dataOut.data_pre[1]
2343 2345
2344 2346 lamb = self.dataOut.C/self.dataOut.frequency
2345 2347 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
2346 2348 paramInterval = self.dataOut.paramInterval
2347 2349
2348 2350 nChannels = data_acf.shape[0]
2349 2351 nLags = data_acf.shape[1]
2350 2352 nProfiles = data_acf.shape[2]
2351 2353 nHeights = self.dataOut.nHeights
2352 2354 nCohInt = self.dataOut.nCohInt
2353 2355 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
2354 2356 heightList = self.dataOut.heightList
2355 2357 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
2356 2358 utctime = self.dataOut.utctime
2357 2359
2358 2360 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
2359 2361
2360 2362 #------------------------ SNR --------------------------------------
2361 2363 power = data_acf[:,0,:,:].real
2362 2364 noise = numpy.zeros(nChannels)
2363 2365 SNR = numpy.zeros(power.shape)
2364 2366 for i in range(nChannels):
2365 2367 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
2366 2368 SNR[i] = (power[i]-noise[i])/noise[i]
2367 2369 SNRm = numpy.nanmean(SNR, axis = 0)
2368 2370 SNRdB = 10*numpy.log10(SNR)
2369 2371
2370 2372 if mode == 'SA':
2371 2373 nPairs = data_ccf.shape[0]
2372 2374 #---------------------- Coherence and Phase --------------------------
2373 2375 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
2374 2376 # phase1 = numpy.copy(phase)
2375 2377 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
2376 2378
2377 2379 for p in range(nPairs):
2378 2380 ch0 = self.dataOut.groupList[p][0]
2379 2381 ch1 = self.dataOut.groupList[p][1]
2380 2382 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
2381 2383 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
2382 2384 # phase1[p,:,:] = numpy.angle(ccf) #median filter
2383 2385 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
2384 2386 # coh1[p,:,:] = numpy.abs(ccf) #median filter
2385 2387 coh = numpy.nanmax(coh1, axis = 0)
2386 2388 # struc = numpy.ones((5,1))
2387 2389 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
2388 2390 #---------------------- Radial Velocity ----------------------------
2389 2391 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
2390 2392 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
2391 2393
2392 2394 if allData:
2393 2395 boolMetFin = ~numpy.isnan(SNRm)
2394 2396 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2395 2397 else:
2396 2398 #------------------------ Meteor mask ---------------------------------
2397 2399 # #SNR mask
2398 2400 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
2399 2401 #
2400 2402 # #Erase small objects
2401 2403 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
2402 2404 #
2403 2405 # auxEEJ = numpy.sum(boolMet1,axis=0)
2404 2406 # indOver = auxEEJ>nProfiles*0.8 #Use this later
2405 2407 # indEEJ = numpy.where(indOver)[0]
2406 2408 # indNEEJ = numpy.where(~indOver)[0]
2407 2409 #
2408 2410 # boolMetFin = boolMet1
2409 2411 #
2410 2412 # if indEEJ.size > 0:
2411 2413 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
2412 2414 #
2413 2415 # boolMet2 = coh > cohThresh
2414 2416 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
2415 2417 #
2416 2418 # #Final Meteor mask
2417 2419 # boolMetFin = boolMet1|boolMet2
2418 2420
2419 2421 #Coherence mask
2420 2422 boolMet1 = coh > 0.75
2421 2423 struc = numpy.ones((30,1))
2422 2424 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
2423 2425
2424 2426 #Derivative mask
2425 2427 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2426 2428 boolMet2 = derPhase < 0.2
2427 2429 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
2428 2430 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
2429 2431 boolMet2 = ndimage.median_filter(boolMet2,size=5)
2430 2432 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
2431 2433 # #Final mask
2432 2434 # boolMetFin = boolMet2
2433 2435 boolMetFin = boolMet1&boolMet2
2434 2436 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
2435 2437 #Creating data_param
2436 2438 coordMet = numpy.where(boolMetFin)
2437 2439
2438 2440 tmet = coordMet[0]
2439 2441 hmet = coordMet[1]
2440 2442
2441 2443 data_param = numpy.zeros((tmet.size, 6 + nPairs))
2442 2444 data_param[:,0] = utctime
2443 2445 data_param[:,1] = tmet
2444 2446 data_param[:,2] = hmet
2445 2447 data_param[:,3] = SNRm[tmet,hmet]
2446 2448 data_param[:,4] = velRad[tmet,hmet]
2447 2449 data_param[:,5] = coh[tmet,hmet]
2448 2450 data_param[:,6:] = phase[:,tmet,hmet].T
2449 2451
2450 2452 elif mode == 'DBS':
2451 2453 self.dataOut.groupList = numpy.arange(nChannels)
2452 2454
2453 2455 #Radial Velocities
2454 2456 # phase = numpy.angle(data_acf[:,1,:,:])
2455 2457 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
2456 2458 velRad = phase*lamb/(4*numpy.pi*tSamp)
2457 2459
2458 2460 #Spectral width
2459 2461 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
2460 2462 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
2461 2463
2462 2464 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
2463 2465 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
2464 2466 if allData:
2465 2467 boolMetFin = ~numpy.isnan(SNRdB)
2466 2468 else:
2467 2469 #SNR
2468 2470 boolMet1 = (SNRdB>SNRthresh) #SNR mask
2469 2471 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
2470 2472
2471 2473 #Radial velocity
2472 2474 boolMet2 = numpy.abs(velRad) < 30
2473 2475 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
2474 2476
2475 2477 #Spectral Width
2476 2478 boolMet3 = spcWidth < 30
2477 2479 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
2478 2480 # boolMetFin = self.__erase_small(boolMet1, 10,5)
2479 2481 boolMetFin = boolMet1&boolMet2&boolMet3
2480 2482
2481 2483 #Creating data_param
2482 2484 coordMet = numpy.where(boolMetFin)
2483 2485
2484 2486 cmet = coordMet[0]
2485 2487 tmet = coordMet[1]
2486 2488 hmet = coordMet[2]
2487 2489
2488 2490 data_param = numpy.zeros((tmet.size, 7))
2489 2491 data_param[:,0] = utctime
2490 2492 data_param[:,1] = cmet
2491 2493 data_param[:,2] = tmet
2492 2494 data_param[:,3] = hmet
2493 2495 data_param[:,4] = SNR[cmet,tmet,hmet].T
2494 2496 data_param[:,5] = velRad[cmet,tmet,hmet].T
2495 2497 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
2496 2498
2497 2499 # self.dataOut.data_param = data_int
2498 2500 if len(data_param) == 0:
2499 2501 self.dataOut.flagNoData = True
2500 2502 else:
2501 2503 self.dataOut.data_param = data_param
2502 2504
2503 2505 def __erase_small(self, binArray, threshX, threshY):
2504 2506 labarray, numfeat = ndimage.measurements.label(binArray)
2505 2507 binArray1 = numpy.copy(binArray)
2506 2508
2507 2509 for i in range(1,numfeat + 1):
2508 2510 auxBin = (labarray==i)
2509 2511 auxSize = auxBin.sum()
2510 2512
2511 2513 x,y = numpy.where(auxBin)
2512 2514 widthX = x.max() - x.min()
2513 2515 widthY = y.max() - y.min()
2514 2516
2515 2517 #width X: 3 seg -> 12.5*3
2516 2518 #width Y:
2517 2519
2518 2520 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
2519 2521 binArray1[auxBin] = False
2520 2522
2521 2523 return binArray1
2522 2524
2523 2525 #--------------- Specular Meteor ----------------
2524 2526
2525 2527 class SMDetection(Operation):
2526 2528 '''
2527 2529 Function DetectMeteors()
2528 2530 Project developed with paper:
2529 2531 HOLDSWORTH ET AL. 2004
2530 2532
2531 2533 Input:
2532 2534 self.dataOut.data_pre
2533 2535
2534 2536 centerReceiverIndex: From the channels, which is the center receiver
2535 2537
2536 2538 hei_ref: Height reference for the Beacon signal extraction
2537 2539 tauindex:
2538 2540 predefinedPhaseShifts: Predefined phase offset for the voltge signals
2539 2541
2540 2542 cohDetection: Whether to user Coherent detection or not
2541 2543 cohDet_timeStep: Coherent Detection calculation time step
2542 2544 cohDet_thresh: Coherent Detection phase threshold to correct phases
2543 2545
2544 2546 noise_timeStep: Noise calculation time step
2545 2547 noise_multiple: Noise multiple to define signal threshold
2546 2548
2547 2549 multDet_timeLimit: Multiple Detection Removal time limit in seconds
2548 2550 multDet_rangeLimit: Multiple Detection Removal range limit in km
2549 2551
2550 2552 phaseThresh: Maximum phase difference between receiver to be consider a meteor
2551 2553 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
2552 2554
2553 2555 hmin: Minimum Height of the meteor to use it in the further wind estimations
2554 2556 hmax: Maximum Height of the meteor to use it in the further wind estimations
2555 2557 azimuth: Azimuth angle correction
2556 2558
2557 2559 Affected:
2558 2560 self.dataOut.data_param
2559 2561
2560 2562 Rejection Criteria (Errors):
2561 2563 0: No error; analysis OK
2562 2564 1: SNR < SNR threshold
2563 2565 2: angle of arrival (AOA) ambiguously determined
2564 2566 3: AOA estimate not feasible
2565 2567 4: Large difference in AOAs obtained from different antenna baselines
2566 2568 5: echo at start or end of time series
2567 2569 6: echo less than 5 examples long; too short for analysis
2568 2570 7: echo rise exceeds 0.3s
2569 2571 8: echo decay time less than twice rise time
2570 2572 9: large power level before echo
2571 2573 10: large power level after echo
2572 2574 11: poor fit to amplitude for estimation of decay time
2573 2575 12: poor fit to CCF phase variation for estimation of radial drift velocity
2574 2576 13: height unresolvable echo: not valid height within 70 to 110 km
2575 2577 14: height ambiguous echo: more then one possible height within 70 to 110 km
2576 2578 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
2577 2579 16: oscilatory echo, indicating event most likely not an underdense echo
2578 2580
2579 2581 17: phase difference in meteor Reestimation
2580 2582
2581 2583 Data Storage:
2582 2584 Meteors for Wind Estimation (8):
2583 2585 Utc Time | Range Height
2584 2586 Azimuth Zenith errorCosDir
2585 2587 VelRad errorVelRad
2586 2588 Phase0 Phase1 Phase2 Phase3
2587 2589 TypeError
2588 2590
2589 2591 '''
2590 2592
2591 2593 def run(self, dataOut, hei_ref = None, tauindex = 0,
2592 2594 phaseOffsets = None,
2593 2595 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
2594 2596 noise_timeStep = 4, noise_multiple = 4,
2595 2597 multDet_timeLimit = 1, multDet_rangeLimit = 3,
2596 2598 phaseThresh = 20, SNRThresh = 5,
2597 2599 hmin = 50, hmax=150, azimuth = 0,
2598 2600 channelPositions = None) :
2599 2601
2600 2602
2601 2603 #Getting Pairslist
2602 2604 if channelPositions == None:
2603 2605 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2604 2606 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2605 2607 meteorOps = SMOperations()
2606 2608 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2607 2609 heiRang = dataOut.getHeiRange()
2608 2610 #Get Beacon signal - No Beacon signal anymore
2609 2611 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
2610 2612 #
2611 2613 # if hei_ref != None:
2612 2614 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
2613 2615 #
2614 2616
2615 2617
2616 2618 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
2617 2619 # see if the user put in pre defined phase shifts
2618 2620 voltsPShift = dataOut.data_pre.copy()
2619 2621
2620 2622 # if predefinedPhaseShifts != None:
2621 2623 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
2622 2624 #
2623 2625 # # elif beaconPhaseShifts:
2624 2626 # # #get hardware phase shifts using beacon signal
2625 2627 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
2626 2628 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
2627 2629 #
2628 2630 # else:
2629 2631 # hardwarePhaseShifts = numpy.zeros(5)
2630 2632 #
2631 2633 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
2632 2634 # for i in range(self.dataOut.data_pre.shape[0]):
2633 2635 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
2634 2636
2635 2637 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
2636 2638
2637 2639 #Remove DC
2638 2640 voltsDC = numpy.mean(voltsPShift,1)
2639 2641 voltsDC = numpy.mean(voltsDC,1)
2640 2642 for i in range(voltsDC.shape[0]):
2641 2643 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
2642 2644
2643 2645 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
2644 2646 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
2645 2647
2646 2648 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
2647 2649 #Coherent Detection
2648 2650 if cohDetection:
2649 2651 #use coherent detection to get the net power
2650 2652 cohDet_thresh = cohDet_thresh*numpy.pi/180
2651 2653 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
2652 2654
2653 2655 #Non-coherent detection!
2654 2656 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
2655 2657 #********** END OF COH/NON-COH POWER CALCULATION**********************
2656 2658
2657 2659 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
2658 2660 #Get noise
2659 2661 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
2660 2662 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
2661 2663 #Get signal threshold
2662 2664 signalThresh = noise_multiple*noise
2663 2665 #Meteor echoes detection
2664 2666 listMeteors = self.__findMeteors(powerNet, signalThresh)
2665 2667 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
2666 2668
2667 2669 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
2668 2670 #Parameters
2669 2671 heiRange = dataOut.getHeiRange()
2670 2672 rangeInterval = heiRange[1] - heiRange[0]
2671 2673 rangeLimit = multDet_rangeLimit/rangeInterval
2672 2674 timeLimit = multDet_timeLimit/dataOut.timeInterval
2673 2675 #Multiple detection removals
2674 2676 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
2675 2677 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
2676 2678
2677 2679 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
2678 2680 #Parameters
2679 2681 phaseThresh = phaseThresh*numpy.pi/180
2680 2682 thresh = [phaseThresh, noise_multiple, SNRThresh]
2681 2683 #Meteor reestimation (Errors N 1, 6, 12, 17)
2682 2684 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
2683 2685 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
2684 2686 #Estimation of decay times (Errors N 7, 8, 11)
2685 2687 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
2686 2688 #******************* END OF METEOR REESTIMATION *******************
2687 2689
2688 2690 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
2689 2691 #Calculating Radial Velocity (Error N 15)
2690 2692 radialStdThresh = 10
2691 2693 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
2692 2694
2693 2695 if len(listMeteors4) > 0:
2694 2696 #Setting New Array
2695 2697 date = dataOut.utctime
2696 2698 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
2697 2699
2698 2700 #Correcting phase offset
2699 2701 if phaseOffsets != None:
2700 2702 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2701 2703 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2702 2704
2703 2705 #Second Pairslist
2704 2706 pairsList = []
2705 2707 pairx = (0,1)
2706 2708 pairy = (2,3)
2707 2709 pairsList.append(pairx)
2708 2710 pairsList.append(pairy)
2709 2711
2710 2712 jph = numpy.array([0,0,0,0])
2711 2713 h = (hmin,hmax)
2712 2714 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2713 2715
2714 2716 # #Calculate AOA (Error N 3, 4)
2715 2717 # #JONES ET AL. 1998
2716 2718 # error = arrayParameters[:,-1]
2717 2719 # AOAthresh = numpy.pi/8
2718 2720 # phases = -arrayParameters[:,9:13]
2719 2721 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
2720 2722 #
2721 2723 # #Calculate Heights (Error N 13 and 14)
2722 2724 # error = arrayParameters[:,-1]
2723 2725 # Ranges = arrayParameters[:,2]
2724 2726 # zenith = arrayParameters[:,5]
2725 2727 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
2726 2728 # error = arrayParameters[:,-1]
2727 2729 #********************* END OF PARAMETERS CALCULATION **************************
2728 2730
2729 2731 #***************************+ PASS DATA TO NEXT STEP **********************
2730 2732 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
2731 2733 dataOut.data_param = arrayParameters
2732 2734
2733 2735 if arrayParameters == None:
2734 2736 dataOut.flagNoData = True
2735 2737 else:
2736 2738 dataOut.flagNoData = True
2737 2739
2738 2740 return
2739 2741
2740 2742 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
2741 2743
2742 2744 minIndex = min(newheis[0])
2743 2745 maxIndex = max(newheis[0])
2744 2746
2745 2747 voltage = voltage0[:,:,minIndex:maxIndex+1]
2746 2748 nLength = voltage.shape[1]/n
2747 2749 nMin = 0
2748 2750 nMax = 0
2749 2751 phaseOffset = numpy.zeros((len(pairslist),n))
2750 2752
2751 2753 for i in range(n):
2752 2754 nMax += nLength
2753 2755 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
2754 2756 phaseCCF = numpy.mean(phaseCCF, axis = 2)
2755 2757 phaseOffset[:,i] = phaseCCF.transpose()
2756 2758 nMin = nMax
2757 2759 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
2758 2760
2759 2761 #Remove Outliers
2760 2762 factor = 2
2761 2763 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
2762 2764 dw = numpy.std(wt,axis = 1)
2763 2765 dw = dw.reshape((dw.size,1))
2764 2766 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
2765 2767 phaseOffset[ind] = numpy.nan
2766 2768 phaseOffset = stats.nanmean(phaseOffset, axis=1)
2767 2769
2768 2770 return phaseOffset
2769 2771
2770 2772 def __shiftPhase(self, data, phaseShift):
2771 2773 #this will shift the phase of a complex number
2772 2774 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
2773 2775 return dataShifted
2774 2776
2775 2777 def __estimatePhaseDifference(self, array, pairslist):
2776 2778 nChannel = array.shape[0]
2777 2779 nHeights = array.shape[2]
2778 2780 numPairs = len(pairslist)
2779 2781 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
2780 2782 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
2781 2783
2782 2784 #Correct phases
2783 2785 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
2784 2786 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2785 2787
2786 2788 if indDer[0].shape[0] > 0:
2787 2789 for i in range(indDer[0].shape[0]):
2788 2790 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
2789 2791 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
2790 2792
2791 2793 # for j in range(numSides):
2792 2794 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
2793 2795 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
2794 2796 #
2795 2797 #Linear
2796 2798 phaseInt = numpy.zeros((numPairs,1))
2797 2799 angAllCCF = phaseCCF[:,[0,1,3,4],0]
2798 2800 for j in range(numPairs):
2799 2801 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
2800 2802 phaseInt[j] = fit[1]
2801 2803 #Phase Differences
2802 2804 phaseDiff = phaseInt - phaseCCF[:,2,:]
2803 2805 phaseArrival = phaseInt.reshape(phaseInt.size)
2804 2806
2805 2807 #Dealias
2806 2808 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
2807 2809 # indAlias = numpy.where(phaseArrival > numpy.pi)
2808 2810 # phaseArrival[indAlias] -= 2*numpy.pi
2809 2811 # indAlias = numpy.where(phaseArrival < -numpy.pi)
2810 2812 # phaseArrival[indAlias] += 2*numpy.pi
2811 2813
2812 2814 return phaseDiff, phaseArrival
2813 2815
2814 2816 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
2815 2817 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
2816 2818 #find the phase shifts of each channel over 1 second intervals
2817 2819 #only look at ranges below the beacon signal
2818 2820 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
2819 2821 numBlocks = int(volts.shape[1]/numProfPerBlock)
2820 2822 numHeights = volts.shape[2]
2821 2823 nChannel = volts.shape[0]
2822 2824 voltsCohDet = volts.copy()
2823 2825
2824 2826 pairsarray = numpy.array(pairslist)
2825 2827 indSides = pairsarray[:,1]
2826 2828 # indSides = numpy.array(range(nChannel))
2827 2829 # indSides = numpy.delete(indSides, indCenter)
2828 2830 #
2829 2831 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
2830 2832 listBlocks = numpy.array_split(volts, numBlocks, 1)
2831 2833
2832 2834 startInd = 0
2833 2835 endInd = 0
2834 2836
2835 2837 for i in range(numBlocks):
2836 2838 startInd = endInd
2837 2839 endInd = endInd + listBlocks[i].shape[1]
2838 2840
2839 2841 arrayBlock = listBlocks[i]
2840 2842 # arrayBlockCenter = listCenter[i]
2841 2843
2842 2844 #Estimate the Phase Difference
2843 2845 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
2844 2846 #Phase Difference RMS
2845 2847 arrayPhaseRMS = numpy.abs(phaseDiff)
2846 2848 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
2847 2849 indPhase = numpy.where(phaseRMSaux==4)
2848 2850 #Shifting
2849 2851 if indPhase[0].shape[0] > 0:
2850 2852 for j in range(indSides.size):
2851 2853 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
2852 2854 voltsCohDet[:,startInd:endInd,:] = arrayBlock
2853 2855
2854 2856 return voltsCohDet
2855 2857
2856 2858 def __calculateCCF(self, volts, pairslist ,laglist):
2857 2859
2858 2860 nHeights = volts.shape[2]
2859 2861 nPoints = volts.shape[1]
2860 2862 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
2861 2863
2862 2864 for i in range(len(pairslist)):
2863 2865 volts1 = volts[pairslist[i][0]]
2864 2866 volts2 = volts[pairslist[i][1]]
2865 2867
2866 2868 for t in range(len(laglist)):
2867 2869 idxT = laglist[t]
2868 2870 if idxT >= 0:
2869 2871 vStacked = numpy.vstack((volts2[idxT:,:],
2870 2872 numpy.zeros((idxT, nHeights),dtype='complex')))
2871 2873 else:
2872 2874 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
2873 2875 volts2[:(nPoints + idxT),:]))
2874 2876 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
2875 2877
2876 2878 vStacked = None
2877 2879 return voltsCCF
2878 2880
2879 2881 def __getNoise(self, power, timeSegment, timeInterval):
2880 2882 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
2881 2883 numBlocks = int(power.shape[0]/numProfPerBlock)
2882 2884 numHeights = power.shape[1]
2883 2885
2884 2886 listPower = numpy.array_split(power, numBlocks, 0)
2885 2887 noise = numpy.zeros((power.shape[0], power.shape[1]))
2886 2888 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
2887 2889
2888 2890 startInd = 0
2889 2891 endInd = 0
2890 2892
2891 2893 for i in range(numBlocks): #split por canal
2892 2894 startInd = endInd
2893 2895 endInd = endInd + listPower[i].shape[0]
2894 2896
2895 2897 arrayBlock = listPower[i]
2896 2898 noiseAux = numpy.mean(arrayBlock, 0)
2897 2899 # noiseAux = numpy.median(noiseAux)
2898 2900 # noiseAux = numpy.mean(arrayBlock)
2899 2901 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
2900 2902
2901 2903 noiseAux1 = numpy.mean(arrayBlock)
2902 2904 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
2903 2905
2904 2906 return noise, noise1
2905 2907
2906 2908 def __findMeteors(self, power, thresh):
2907 2909 nProf = power.shape[0]
2908 2910 nHeights = power.shape[1]
2909 2911 listMeteors = []
2910 2912
2911 2913 for i in range(nHeights):
2912 2914 powerAux = power[:,i]
2913 2915 threshAux = thresh[:,i]
2914 2916
2915 2917 indUPthresh = numpy.where(powerAux > threshAux)[0]
2916 2918 indDNthresh = numpy.where(powerAux <= threshAux)[0]
2917 2919
2918 2920 j = 0
2919 2921
2920 2922 while (j < indUPthresh.size - 2):
2921 2923 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
2922 2924 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
2923 2925 indDNthresh = indDNthresh[indDNAux]
2924 2926
2925 2927 if (indDNthresh.size > 0):
2926 2928 indEnd = indDNthresh[0] - 1
2927 2929 indInit = indUPthresh[j]
2928 2930
2929 2931 meteor = powerAux[indInit:indEnd + 1]
2930 2932 indPeak = meteor.argmax() + indInit
2931 2933 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
2932 2934
2933 2935 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
2934 2936 j = numpy.where(indUPthresh == indEnd)[0] + 1
2935 2937 else: j+=1
2936 2938 else: j+=1
2937 2939
2938 2940 return listMeteors
2939 2941
2940 2942 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
2941 2943
2942 2944 arrayMeteors = numpy.asarray(listMeteors)
2943 2945 listMeteors1 = []
2944 2946
2945 2947 while arrayMeteors.shape[0] > 0:
2946 2948 FLAs = arrayMeteors[:,4]
2947 2949 maxFLA = FLAs.argmax()
2948 2950 listMeteors1.append(arrayMeteors[maxFLA,:])
2949 2951
2950 2952 MeteorInitTime = arrayMeteors[maxFLA,1]
2951 2953 MeteorEndTime = arrayMeteors[maxFLA,3]
2952 2954 MeteorHeight = arrayMeteors[maxFLA,0]
2953 2955
2954 2956 #Check neighborhood
2955 2957 maxHeightIndex = MeteorHeight + rangeLimit
2956 2958 minHeightIndex = MeteorHeight - rangeLimit
2957 2959 minTimeIndex = MeteorInitTime - timeLimit
2958 2960 maxTimeIndex = MeteorEndTime + timeLimit
2959 2961
2960 2962 #Check Heights
2961 2963 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
2962 2964 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
2963 2965 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
2964 2966
2965 2967 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
2966 2968
2967 2969 return listMeteors1
2968 2970
2969 2971 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
2970 2972 numHeights = volts.shape[2]
2971 2973 nChannel = volts.shape[0]
2972 2974
2973 2975 thresholdPhase = thresh[0]
2974 2976 thresholdNoise = thresh[1]
2975 2977 thresholdDB = float(thresh[2])
2976 2978
2977 2979 thresholdDB1 = 10**(thresholdDB/10)
2978 2980 pairsarray = numpy.array(pairslist)
2979 2981 indSides = pairsarray[:,1]
2980 2982
2981 2983 pairslist1 = list(pairslist)
2982 2984 pairslist1.append((0,1))
2983 2985 pairslist1.append((3,4))
2984 2986
2985 2987 listMeteors1 = []
2986 2988 listPowerSeries = []
2987 2989 listVoltageSeries = []
2988 2990 #volts has the war data
2989 2991
2990 2992 if frequency == 30e6:
2991 2993 timeLag = 45*10**-3
2992 2994 else:
2993 2995 timeLag = 15*10**-3
2994 2996 lag = numpy.ceil(timeLag/timeInterval)
2995 2997
2996 2998 for i in range(len(listMeteors)):
2997 2999
2998 3000 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
2999 3001 meteorAux = numpy.zeros(16)
3000 3002
3001 3003 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
3002 3004 mHeight = listMeteors[i][0]
3003 3005 mStart = listMeteors[i][1]
3004 3006 mPeak = listMeteors[i][2]
3005 3007 mEnd = listMeteors[i][3]
3006 3008
3007 3009 #get the volt data between the start and end times of the meteor
3008 3010 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
3009 3011 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3010 3012
3011 3013 #3.6. Phase Difference estimation
3012 3014 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
3013 3015
3014 3016 #3.7. Phase difference removal & meteor start, peak and end times reestimated
3015 3017 #meteorVolts0.- all Channels, all Profiles
3016 3018 meteorVolts0 = volts[:,:,mHeight]
3017 3019 meteorThresh = noise[:,mHeight]*thresholdNoise
3018 3020 meteorNoise = noise[:,mHeight]
3019 3021 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
3020 3022 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
3021 3023
3022 3024 #Times reestimation
3023 3025 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
3024 3026 if mStart1.size > 0:
3025 3027 mStart1 = mStart1[-1] + 1
3026 3028
3027 3029 else:
3028 3030 mStart1 = mPeak
3029 3031
3030 3032 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
3031 3033 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
3032 3034 if mEndDecayTime1.size == 0:
3033 3035 mEndDecayTime1 = powerNet0.size
3034 3036 else:
3035 3037 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
3036 3038 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
3037 3039
3038 3040 #meteorVolts1.- all Channels, from start to end
3039 3041 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
3040 3042 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
3041 3043 if meteorVolts2.shape[1] == 0:
3042 3044 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
3043 3045 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
3044 3046 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
3045 3047 ##################### END PARAMETERS REESTIMATION #########################
3046 3048
3047 3049 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
3048 3050 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
3049 3051 if meteorVolts2.shape[1] > 0:
3050 3052 #Phase Difference re-estimation
3051 3053 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
3052 3054 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
3053 3055 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
3054 3056 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
3055 3057 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
3056 3058
3057 3059 #Phase Difference RMS
3058 3060 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
3059 3061 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
3060 3062 #Data from Meteor
3061 3063 mPeak1 = powerNet1.argmax() + mStart1
3062 3064 mPeakPower1 = powerNet1.max()
3063 3065 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
3064 3066 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
3065 3067 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
3066 3068 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
3067 3069 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
3068 3070 #Vectorize
3069 3071 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
3070 3072 meteorAux[7:11] = phaseDiffint[0:4]
3071 3073
3072 3074 #Rejection Criterions
3073 3075 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
3074 3076 meteorAux[-1] = 17
3075 3077 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
3076 3078 meteorAux[-1] = 1
3077 3079
3078 3080
3079 3081 else:
3080 3082 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
3081 3083 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
3082 3084 PowerSeries = 0
3083 3085
3084 3086 listMeteors1.append(meteorAux)
3085 3087 listPowerSeries.append(PowerSeries)
3086 3088 listVoltageSeries.append(meteorVolts1)
3087 3089
3088 3090 return listMeteors1, listPowerSeries, listVoltageSeries
3089 3091
3090 3092 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
3091 3093
3092 3094 threshError = 10
3093 3095 #Depending if it is 30 or 50 MHz
3094 3096 if frequency == 30e6:
3095 3097 timeLag = 45*10**-3
3096 3098 else:
3097 3099 timeLag = 15*10**-3
3098 3100 lag = numpy.ceil(timeLag/timeInterval)
3099 3101
3100 3102 listMeteors1 = []
3101 3103
3102 3104 for i in range(len(listMeteors)):
3103 3105 meteorPower = listPower[i]
3104 3106 meteorAux = listMeteors[i]
3105 3107
3106 3108 if meteorAux[-1] == 0:
3107 3109
3108 3110 try:
3109 3111 indmax = meteorPower.argmax()
3110 3112 indlag = indmax + lag
3111 3113
3112 3114 y = meteorPower[indlag:]
3113 3115 x = numpy.arange(0, y.size)*timeLag
3114 3116
3115 3117 #first guess
3116 3118 a = y[0]
3117 3119 tau = timeLag
3118 3120 #exponential fit
3119 3121 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
3120 3122 y1 = self.__exponential_function(x, *popt)
3121 3123 #error estimation
3122 3124 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
3123 3125
3124 3126 decayTime = popt[1]
3125 3127 riseTime = indmax*timeInterval
3126 3128 meteorAux[11:13] = [decayTime, error]
3127 3129
3128 3130 #Table items 7, 8 and 11
3129 3131 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
3130 3132 meteorAux[-1] = 7
3131 3133 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
3132 3134 meteorAux[-1] = 8
3133 3135 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
3134 3136 meteorAux[-1] = 11
3135 3137
3136 3138
3137 3139 except:
3138 3140 meteorAux[-1] = 11
3139 3141
3140 3142
3141 3143 listMeteors1.append(meteorAux)
3142 3144
3143 3145 return listMeteors1
3144 3146
3145 3147 #Exponential Function
3146 3148
3147 3149 def __exponential_function(self, x, a, tau):
3148 3150 y = a*numpy.exp(-x/tau)
3149 3151 return y
3150 3152
3151 3153 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
3152 3154
3153 3155 pairslist1 = list(pairslist)
3154 3156 pairslist1.append((0,1))
3155 3157 pairslist1.append((3,4))
3156 3158 numPairs = len(pairslist1)
3157 3159 #Time Lag
3158 3160 timeLag = 45*10**-3
3159 3161 c = 3e8
3160 3162 lag = numpy.ceil(timeLag/timeInterval)
3161 3163 freq = 30e6
3162 3164
3163 3165 listMeteors1 = []
3164 3166
3165 3167 for i in range(len(listMeteors)):
3166 3168 meteorAux = listMeteors[i]
3167 3169 if meteorAux[-1] == 0:
3168 3170 mStart = listMeteors[i][1]
3169 3171 mPeak = listMeteors[i][2]
3170 3172 mLag = mPeak - mStart + lag
3171 3173
3172 3174 #get the volt data between the start and end times of the meteor
3173 3175 meteorVolts = listVolts[i]
3174 3176 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3175 3177
3176 3178 #Get CCF
3177 3179 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
3178 3180
3179 3181 #Method 2
3180 3182 slopes = numpy.zeros(numPairs)
3181 3183 time = numpy.array([-2,-1,1,2])*timeInterval
3182 3184 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
3183 3185
3184 3186 #Correct phases
3185 3187 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
3186 3188 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
3187 3189
3188 3190 if indDer[0].shape[0] > 0:
3189 3191 for i in range(indDer[0].shape[0]):
3190 3192 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
3191 3193 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
3192 3194
3193 3195 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
3194 3196 for j in range(numPairs):
3195 3197 fit = stats.linregress(time, angAllCCF[j,:])
3196 3198 slopes[j] = fit[0]
3197 3199
3198 3200 #Remove Outlier
3199 3201 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3200 3202 # slopes = numpy.delete(slopes,indOut)
3201 3203 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3202 3204 # slopes = numpy.delete(slopes,indOut)
3203 3205
3204 3206 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
3205 3207 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
3206 3208 meteorAux[-2] = radialError
3207 3209 meteorAux[-3] = radialVelocity
3208 3210
3209 3211 #Setting Error
3210 3212 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
3211 3213 if numpy.abs(radialVelocity) > 200:
3212 3214 meteorAux[-1] = 15
3213 3215 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
3214 3216 elif radialError > radialStdThresh:
3215 3217 meteorAux[-1] = 12
3216 3218
3217 3219 listMeteors1.append(meteorAux)
3218 3220 return listMeteors1
3219 3221
3220 3222 def __setNewArrays(self, listMeteors, date, heiRang):
3221 3223
3222 3224 #New arrays
3223 3225 arrayMeteors = numpy.array(listMeteors)
3224 3226 arrayParameters = numpy.zeros((len(listMeteors), 13))
3225 3227
3226 3228 #Date inclusion
3227 3229 # date = re.findall(r'\((.*?)\)', date)
3228 3230 # date = date[0].split(',')
3229 3231 # date = map(int, date)
3230 3232 #
3231 3233 # if len(date)<6:
3232 3234 # date.append(0)
3233 3235 #
3234 3236 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
3235 3237 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
3236 3238 arrayDate = numpy.tile(date, (len(listMeteors)))
3237 3239
3238 3240 #Meteor array
3239 3241 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
3240 3242 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
3241 3243
3242 3244 #Parameters Array
3243 3245 arrayParameters[:,0] = arrayDate #Date
3244 3246 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
3245 3247 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
3246 3248 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
3247 3249 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
3248 3250
3249 3251
3250 3252 return arrayParameters
3251 3253
3252 3254 class CorrectSMPhases(Operation):
3253 3255
3254 3256 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
3255 3257
3256 3258 arrayParameters = dataOut.data_param
3257 3259 pairsList = []
3258 3260 pairx = (0,1)
3259 3261 pairy = (2,3)
3260 3262 pairsList.append(pairx)
3261 3263 pairsList.append(pairy)
3262 3264 jph = numpy.zeros(4)
3263 3265
3264 3266 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
3265 3267 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
3266 3268 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
3267 3269
3268 3270 meteorOps = SMOperations()
3269 3271 if channelPositions == None:
3270 3272 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3271 3273 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3272 3274
3273 3275 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3274 3276 h = (hmin,hmax)
3275 3277
3276 3278 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
3277 3279
3278 3280 dataOut.data_param = arrayParameters
3279 3281 return
3280 3282
3281 3283 class SMPhaseCalibration(Operation):
3282 3284
3283 3285 __buffer = None
3284 3286
3285 3287 __initime = None
3286 3288
3287 3289 __dataReady = False
3288 3290
3289 3291 __isConfig = False
3290 3292
3291 3293 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
3292 3294
3293 3295 dataTime = currentTime + paramInterval
3294 3296 deltaTime = dataTime - initTime
3295 3297
3296 3298 if deltaTime >= outputInterval or deltaTime < 0:
3297 3299 return True
3298 3300
3299 3301 return False
3300 3302
3301 3303 def __getGammas(self, pairs, d, phases):
3302 3304 gammas = numpy.zeros(2)
3303 3305
3304 3306 for i in range(len(pairs)):
3305 3307
3306 3308 pairi = pairs[i]
3307 3309
3308 3310 phip3 = phases[:,pairi[1]]
3309 3311 d3 = d[pairi[1]]
3310 3312 phip2 = phases[:,pairi[0]]
3311 3313 d2 = d[pairi[0]]
3312 3314 #Calculating gamma
3313 3315 # jdcos = alp1/(k*d1)
3314 3316 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
3315 3317 jgamma = -phip2*d3/d2 - phip3
3316 3318 jgamma = numpy.angle(numpy.exp(1j*jgamma))
3317 3319 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
3318 3320 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
3319 3321
3320 3322 #Revised distribution
3321 3323 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
3322 3324
3323 3325 #Histogram
3324 3326 nBins = 64.0
3325 3327 rmin = -0.5*numpy.pi
3326 3328 rmax = 0.5*numpy.pi
3327 3329 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
3328 3330
3329 3331 meteorsY = phaseHisto[0]
3330 3332 phasesX = phaseHisto[1][:-1]
3331 3333 width = phasesX[1] - phasesX[0]
3332 3334 phasesX += width/2
3333 3335
3334 3336 #Gaussian aproximation
3335 3337 bpeak = meteorsY.argmax()
3336 3338 peak = meteorsY.max()
3337 3339 jmin = bpeak - 5
3338 3340 jmax = bpeak + 5 + 1
3339 3341
3340 3342 if jmin<0:
3341 3343 jmin = 0
3342 3344 jmax = 6
3343 3345 elif jmax > meteorsY.size:
3344 3346 jmin = meteorsY.size - 6
3345 3347 jmax = meteorsY.size
3346 3348
3347 3349 x0 = numpy.array([peak,bpeak,50])
3348 3350 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
3349 3351
3350 3352 #Gammas
3351 3353 gammas[i] = coeff[0][1]
3352 3354
3353 3355 return gammas
3354 3356
3355 3357 def __residualFunction(self, coeffs, y, t):
3356 3358
3357 3359 return y - self.__gauss_function(t, coeffs)
3358 3360
3359 3361 def __gauss_function(self, t, coeffs):
3360 3362
3361 3363 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
3362 3364
3363 3365 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
3364 3366 meteorOps = SMOperations()
3365 3367 nchan = 4
3366 3368 pairx = pairsList[0]
3367 3369 pairy = pairsList[1]
3368 3370 center_xangle = 0
3369 3371 center_yangle = 0
3370 3372 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
3371 3373 ntimes = len(range_angle)
3372 3374
3373 3375 nstepsx = 20.0
3374 3376 nstepsy = 20.0
3375 3377
3376 3378 for iz in range(ntimes):
3377 3379 min_xangle = -range_angle[iz]/2 + center_xangle
3378 3380 max_xangle = range_angle[iz]/2 + center_xangle
3379 3381 min_yangle = -range_angle[iz]/2 + center_yangle
3380 3382 max_yangle = range_angle[iz]/2 + center_yangle
3381 3383
3382 3384 inc_x = (max_xangle-min_xangle)/nstepsx
3383 3385 inc_y = (max_yangle-min_yangle)/nstepsy
3384 3386
3385 3387 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
3386 3388 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
3387 3389 penalty = numpy.zeros((nstepsx,nstepsy))
3388 3390 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
3389 3391 jph = numpy.zeros(nchan)
3390 3392
3391 3393 # Iterations looking for the offset
3392 3394 for iy in range(int(nstepsy)):
3393 3395 for ix in range(int(nstepsx)):
3394 3396 jph[pairy[1]] = alpha_y[iy]
3395 3397 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
3396 3398
3397 3399 jph[pairx[1]] = alpha_x[ix]
3398 3400 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
3399 3401
3400 3402 jph_array[:,ix,iy] = jph
3401 3403
3402 3404 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
3403 3405 error = meteorsArray1[:,-1]
3404 3406 ind1 = numpy.where(error==0)[0]
3405 3407 penalty[ix,iy] = ind1.size
3406 3408
3407 3409 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
3408 3410 phOffset = jph_array[:,i,j]
3409 3411
3410 3412 center_xangle = phOffset[pairx[1]]
3411 3413 center_yangle = phOffset[pairy[1]]
3412 3414
3413 3415 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
3414 3416 phOffset = phOffset*180/numpy.pi
3415 3417 return phOffset
3416 3418
3417 3419
3418 3420 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
3419 3421
3420 3422 dataOut.flagNoData = True
3421 3423 self.__dataReady = False
3422 3424 dataOut.outputInterval = nHours*3600
3423 3425
3424 3426 if self.__isConfig == False:
3425 3427 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
3426 3428 #Get Initial LTC time
3427 3429 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3428 3430 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
3429 3431
3430 3432 self.__isConfig = True
3431 3433
3432 3434 if self.__buffer == None:
3433 3435 self.__buffer = dataOut.data_param.copy()
3434 3436
3435 3437 else:
3436 3438 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
3437 3439
3438 3440 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
3439 3441
3440 3442 if self.__dataReady:
3441 3443 dataOut.utctimeInit = self.__initime
3442 3444 self.__initime += dataOut.outputInterval #to erase time offset
3443 3445
3444 3446 freq = dataOut.frequency
3445 3447 c = dataOut.C #m/s
3446 3448 lamb = c/freq
3447 3449 k = 2*numpy.pi/lamb
3448 3450 azimuth = 0
3449 3451 h = (hmin, hmax)
3450 3452 pairs = ((0,1),(2,3))
3451 3453
3452 3454 if channelPositions == None:
3453 3455 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3454 3456 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3455 3457 meteorOps = SMOperations()
3456 3458 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3457 3459
3458 3460 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
3459 3461
3460 3462 meteorsArray = self.__buffer
3461 3463 error = meteorsArray[:,-1]
3462 3464 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
3463 3465 ind1 = numpy.where(boolError)[0]
3464 3466 meteorsArray = meteorsArray[ind1,:]
3465 3467 meteorsArray[:,-1] = 0
3466 3468 phases = meteorsArray[:,8:12]
3467 3469
3468 3470 #Calculate Gammas
3469 3471 gammas = self.__getGammas(pairs, distances, phases)
3470 3472 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
3471 3473 #Calculate Phases
3472 3474 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
3473 3475 phasesOff = phasesOff.reshape((1,phasesOff.size))
3474 3476 dataOut.data_output = -phasesOff
3475 3477 dataOut.flagNoData = False
3476 3478 self.__buffer = None
3477 3479
3478 3480
3479 3481 return
3480 3482
3481 3483 class SMOperations():
3482 3484
3483 3485 def __init__(self):
3484 3486
3485 3487 return
3486 3488
3487 3489 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
3488 3490
3489 3491 arrayParameters = arrayParameters0.copy()
3490 3492 hmin = h[0]
3491 3493 hmax = h[1]
3492 3494
3493 3495 #Calculate AOA (Error N 3, 4)
3494 3496 #JONES ET AL. 1998
3495 3497 AOAthresh = numpy.pi/8
3496 3498 error = arrayParameters[:,-1]
3497 3499 phases = -arrayParameters[:,8:12] + jph
3498 3500 # phases = numpy.unwrap(phases)
3499 3501 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
3500 3502
3501 3503 #Calculate Heights (Error N 13 and 14)
3502 3504 error = arrayParameters[:,-1]
3503 3505 Ranges = arrayParameters[:,1]
3504 3506 zenith = arrayParameters[:,4]
3505 3507 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
3506 3508
3507 3509 #----------------------- Get Final data ------------------------------------
3508 3510 # error = arrayParameters[:,-1]
3509 3511 # ind1 = numpy.where(error==0)[0]
3510 3512 # arrayParameters = arrayParameters[ind1,:]
3511 3513
3512 3514 return arrayParameters
3513 3515
3514 3516 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
3515 3517
3516 3518 arrayAOA = numpy.zeros((phases.shape[0],3))
3517 3519 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
3518 3520
3519 3521 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3520 3522 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3521 3523 arrayAOA[:,2] = cosDirError
3522 3524
3523 3525 azimuthAngle = arrayAOA[:,0]
3524 3526 zenithAngle = arrayAOA[:,1]
3525 3527
3526 3528 #Setting Error
3527 3529 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
3528 3530 error[indError] = 0
3529 3531 #Number 3: AOA not fesible
3530 3532 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3531 3533 error[indInvalid] = 3
3532 3534 #Number 4: Large difference in AOAs obtained from different antenna baselines
3533 3535 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3534 3536 error[indInvalid] = 4
3535 3537 return arrayAOA, error
3536 3538
3537 3539 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
3538 3540
3539 3541 #Initializing some variables
3540 3542 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3541 3543 ang_aux = ang_aux.reshape(1,ang_aux.size)
3542 3544
3543 3545 cosdir = numpy.zeros((arrayPhase.shape[0],2))
3544 3546 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3545 3547
3546 3548
3547 3549 for i in range(2):
3548 3550 ph0 = arrayPhase[:,pairsList[i][0]]
3549 3551 ph1 = arrayPhase[:,pairsList[i][1]]
3550 3552 d0 = distances[pairsList[i][0]]
3551 3553 d1 = distances[pairsList[i][1]]
3552 3554
3553 3555 ph0_aux = ph0 + ph1
3554 3556 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
3555 3557 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
3556 3558 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
3557 3559 #First Estimation
3558 3560 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
3559 3561
3560 3562 #Most-Accurate Second Estimation
3561 3563 phi1_aux = ph0 - ph1
3562 3564 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3563 3565 #Direction Cosine 1
3564 3566 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
3565 3567
3566 3568 #Searching the correct Direction Cosine
3567 3569 cosdir0_aux = cosdir0[:,i]
3568 3570 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3569 3571 #Minimum Distance
3570 3572 cosDiff = (cosdir1 - cosdir0_aux)**2
3571 3573 indcos = cosDiff.argmin(axis = 1)
3572 3574 #Saving Value obtained
3573 3575 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3574 3576
3575 3577 return cosdir0, cosdir
3576 3578
3577 3579 def __calculateAOA(self, cosdir, azimuth):
3578 3580 cosdirX = cosdir[:,0]
3579 3581 cosdirY = cosdir[:,1]
3580 3582
3581 3583 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3582 3584 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
3583 3585 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3584 3586
3585 3587 return angles
3586 3588
3587 3589 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3588 3590
3589 3591 Ramb = 375 #Ramb = c/(2*PRF)
3590 3592 Re = 6371 #Earth Radius
3591 3593 heights = numpy.zeros(Ranges.shape)
3592 3594
3593 3595 R_aux = numpy.array([0,1,2])*Ramb
3594 3596 R_aux = R_aux.reshape(1,R_aux.size)
3595 3597
3596 3598 Ranges = Ranges.reshape(Ranges.size,1)
3597 3599
3598 3600 Ri = Ranges + R_aux
3599 3601 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3600 3602
3601 3603 #Check if there is a height between 70 and 110 km
3602 3604 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3603 3605 ind_h = numpy.where(h_bool == 1)[0]
3604 3606
3605 3607 hCorr = hi[ind_h, :]
3606 3608 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3607 3609
3608 3610 hCorr = hi[ind_hCorr]
3609 3611 heights[ind_h] = hCorr
3610 3612
3611 3613 #Setting Error
3612 3614 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3613 3615 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3614 3616 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
3615 3617 error[indError] = 0
3616 3618 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3617 3619 error[indInvalid2] = 14
3618 3620 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3619 3621 error[indInvalid1] = 13
3620 3622
3621 3623 return heights, error
3622 3624
3623 3625 def getPhasePairs(self, channelPositions):
3624 3626 chanPos = numpy.array(channelPositions)
3625 3627 listOper = list(itertools.combinations(range(5),2))
3626 3628
3627 3629 distances = numpy.zeros(4)
3628 3630 axisX = []
3629 3631 axisY = []
3630 3632 distX = numpy.zeros(3)
3631 3633 distY = numpy.zeros(3)
3632 3634 ix = 0
3633 3635 iy = 0
3634 3636
3635 3637 pairX = numpy.zeros((2,2))
3636 3638 pairY = numpy.zeros((2,2))
3637 3639
3638 3640 for i in range(len(listOper)):
3639 3641 pairi = listOper[i]
3640 3642
3641 3643 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
3642 3644
3643 3645 if posDif[0] == 0:
3644 3646 axisY.append(pairi)
3645 3647 distY[iy] = posDif[1]
3646 3648 iy += 1
3647 3649 elif posDif[1] == 0:
3648 3650 axisX.append(pairi)
3649 3651 distX[ix] = posDif[0]
3650 3652 ix += 1
3651 3653
3652 3654 for i in range(2):
3653 3655 if i==0:
3654 3656 dist0 = distX
3655 3657 axis0 = axisX
3656 3658 else:
3657 3659 dist0 = distY
3658 3660 axis0 = axisY
3659 3661
3660 3662 side = numpy.argsort(dist0)[:-1]
3661 3663 axis0 = numpy.array(axis0)[side,:]
3662 3664 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
3663 3665 axis1 = numpy.unique(numpy.reshape(axis0,4))
3664 3666 side = axis1[axis1 != chanC]
3665 3667 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
3666 3668 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
3667 3669 if diff1<0:
3668 3670 chan2 = side[0]
3669 3671 d2 = numpy.abs(diff1)
3670 3672 chan1 = side[1]
3671 3673 d1 = numpy.abs(diff2)
3672 3674 else:
3673 3675 chan2 = side[1]
3674 3676 d2 = numpy.abs(diff2)
3675 3677 chan1 = side[0]
3676 3678 d1 = numpy.abs(diff1)
3677 3679
3678 3680 if i==0:
3679 3681 chanCX = chanC
3680 3682 chan1X = chan1
3681 3683 chan2X = chan2
3682 3684 distances[0:2] = numpy.array([d1,d2])
3683 3685 else:
3684 3686 chanCY = chanC
3685 3687 chan1Y = chan1
3686 3688 chan2Y = chan2
3687 3689 distances[2:4] = numpy.array([d1,d2])
3688 3690 # axisXsides = numpy.reshape(axisX[ix,:],4)
3689 3691 #
3690 3692 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
3691 3693 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
3692 3694 #
3693 3695 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
3694 3696 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
3695 3697 # channel25X = int(pairX[0,ind25X])
3696 3698 # channel20X = int(pairX[1,ind20X])
3697 3699 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
3698 3700 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
3699 3701 # channel25Y = int(pairY[0,ind25Y])
3700 3702 # channel20Y = int(pairY[1,ind20Y])
3701 3703
3702 3704 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
3703 3705 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
3704 3706
3705 3707 return pairslist, distances
3706 3708 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
3707 3709 #
3708 3710 # arrayAOA = numpy.zeros((phases.shape[0],3))
3709 3711 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
3710 3712 #
3711 3713 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3712 3714 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3713 3715 # arrayAOA[:,2] = cosDirError
3714 3716 #
3715 3717 # azimuthAngle = arrayAOA[:,0]
3716 3718 # zenithAngle = arrayAOA[:,1]
3717 3719 #
3718 3720 # #Setting Error
3719 3721 # #Number 3: AOA not fesible
3720 3722 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3721 3723 # error[indInvalid] = 3
3722 3724 # #Number 4: Large difference in AOAs obtained from different antenna baselines
3723 3725 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3724 3726 # error[indInvalid] = 4
3725 3727 # return arrayAOA, error
3726 3728 #
3727 3729 # def __getDirectionCosines(self, arrayPhase, pairsList):
3728 3730 #
3729 3731 # #Initializing some variables
3730 3732 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3731 3733 # ang_aux = ang_aux.reshape(1,ang_aux.size)
3732 3734 #
3733 3735 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
3734 3736 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3735 3737 #
3736 3738 #
3737 3739 # for i in range(2):
3738 3740 # #First Estimation
3739 3741 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
3740 3742 # #Dealias
3741 3743 # indcsi = numpy.where(phi0_aux > numpy.pi)
3742 3744 # phi0_aux[indcsi] -= 2*numpy.pi
3743 3745 # indcsi = numpy.where(phi0_aux < -numpy.pi)
3744 3746 # phi0_aux[indcsi] += 2*numpy.pi
3745 3747 # #Direction Cosine 0
3746 3748 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
3747 3749 #
3748 3750 # #Most-Accurate Second Estimation
3749 3751 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
3750 3752 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3751 3753 # #Direction Cosine 1
3752 3754 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
3753 3755 #
3754 3756 # #Searching the correct Direction Cosine
3755 3757 # cosdir0_aux = cosdir0[:,i]
3756 3758 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3757 3759 # #Minimum Distance
3758 3760 # cosDiff = (cosdir1 - cosdir0_aux)**2
3759 3761 # indcos = cosDiff.argmin(axis = 1)
3760 3762 # #Saving Value obtained
3761 3763 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3762 3764 #
3763 3765 # return cosdir0, cosdir
3764 3766 #
3765 3767 # def __calculateAOA(self, cosdir, azimuth):
3766 3768 # cosdirX = cosdir[:,0]
3767 3769 # cosdirY = cosdir[:,1]
3768 3770 #
3769 3771 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3770 3772 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
3771 3773 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3772 3774 #
3773 3775 # return angles
3774 3776 #
3775 3777 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3776 3778 #
3777 3779 # Ramb = 375 #Ramb = c/(2*PRF)
3778 3780 # Re = 6371 #Earth Radius
3779 3781 # heights = numpy.zeros(Ranges.shape)
3780 3782 #
3781 3783 # R_aux = numpy.array([0,1,2])*Ramb
3782 3784 # R_aux = R_aux.reshape(1,R_aux.size)
3783 3785 #
3784 3786 # Ranges = Ranges.reshape(Ranges.size,1)
3785 3787 #
3786 3788 # Ri = Ranges + R_aux
3787 3789 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3788 3790 #
3789 3791 # #Check if there is a height between 70 and 110 km
3790 3792 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3791 3793 # ind_h = numpy.where(h_bool == 1)[0]
3792 3794 #
3793 3795 # hCorr = hi[ind_h, :]
3794 3796 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3795 3797 #
3796 3798 # hCorr = hi[ind_hCorr]
3797 3799 # heights[ind_h] = hCorr
3798 3800 #
3799 3801 # #Setting Error
3800 3802 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3801 3803 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3802 3804 #
3803 3805 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3804 3806 # error[indInvalid2] = 14
3805 3807 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3806 3808 # error[indInvalid1] = 13
3807 3809 #
3808 3810 # return heights, error
3809 3811 No newline at end of file
@@ -1,1 +1,1
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /><Parameter format="date" id="191113" name="startDate" value="2018/02/23" /><Parameter format="date" id="191114" name="endDate" value="2018/02/23" /><Parameter format="time" id="191115" name="startTime" value="00:00:00" /><Parameter format="time" id="191116" name="endTime" value="03:59:50" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="PositiveLimit" value="1.5" /><Parameter format="float" id="191322" name="NegativeLimit" value="12.5" /></Operation><Operation id="19133" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191331" name="SNRlimit" value="-16" /><Parameter format="float" id="191332" name="Xi01" value="1.500" /><Parameter format="float" id="191333" name="Xi02" value="1.500" /><Parameter format="float" id="191334" name="Xi12" value="0" /><Parameter format="float" id="191335" name="Eta01" value="0.875" /><Parameter format="float" id="191336" name="Eta02" value="-0.875" /><Parameter format="float" id="191337" name="Eta12" value="-1.750" /></Operation><Operation id="19134" name="WindProfilerPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="4" /><Parameter format="str" id="191342" name="wintitle" value="Wind Profiler" /><Parameter format="float" id="191343" name="xmin" value="0" /><Parameter format="float" id="191344" name="xmax" value="4" /><Parameter format="float" id="191345" name="ymin" value="0" /><Parameter format="int" id="191346" name="ymax" value="4" /><Parameter format="float" id="191347" name="zmin" value="-20" /><Parameter format="float" id="191348" name="zmax" value="20" /><Parameter format="float" id="191349" name="SNRmin" value="-20" /><Parameter format="float" id="191350" name="SNRmax" value="20" /><Parameter format="float" id="191351" name="zmin_ver" value="-200" /><Parameter format="float" id="191352" name="zmax_ver" value="200" /><Parameter format="float" id="191353" name="SNRthresh" value="-20" /><Parameter format="int" id="191354" name="save" value="1" /><Parameter format="str" id="191355" name="figpath" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/ImagenesTesis" /></Operation><Operation id="19135" name="ParamWriter" priority="5" type="other"><Parameter format="str" id="191351" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/ImagenesTesis" /><Parameter format="int" id="191352" name="blocksPerFile" value="500" /><Parameter format="list" id="191353" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191354" name="dataList" value="data_output,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation><Operation id="19123" name="setH0" priority="3" type="self"><Parameter format="float" id="191231" name="h0" value="0.500" /></Operation></ProcUnit></Project> No newline at end of file
1 <Project description="ProcBLTR Test" id="191" name="test01"><ReadUnit datatype="BLTRSpectraReader" id="1911" inputId="0" name="BLTRSpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="BLTRSpectraReader" /><Parameter format="str" id="191112" name="path" value="/data/BLTR/2018/enero" /><Parameter format="date" id="191113" name="startDate" value="2016/11/01" /><Parameter format="date" id="191114" name="endDate" value="2016/11/01" /><Parameter format="time" id="191115" name="startTime" value="0:00:00" /><Parameter format="time" id="191116" name="endTime" value="23:59:59" /><Parameter format="int" id="191118" name="walk" value="0" /><Parameter format="str" id="191119" name="ReadMode" value="1" /><Parameter format="int" id="191120" name="online" value="0" /></Operation></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="FullSpectralAnalysis" priority="2" type="other"><Parameter format="float" id="191321" name="SNRlimit" value="2" /></Operation><Operation id="19133" name="WindProfilerPlot" priority="3" type="other"><Parameter format="int" id="191331" name="id" value="4" /><Parameter format="str" id="191332" name="wintitle" value="Wind Profiler" /><Parameter format="int" id="191333" name="zmin" value="-20" /><Parameter format="int" id="191334" name="zmax" value="20" /><Parameter format="float" id="191335" name="zmin_ver" value="-300" /><Parameter format="float" id="191336" name="zmax_ver" value="300" /><Parameter format="int" id="191337" name="SNRmin" value="-5" /><Parameter format="int" id="191338" name="SNRmax" value="30" /><Parameter format="float" id="191339" name="xmin" value="4" /><Parameter format="float" id="191340" name="xmax" value="9" /><Parameter format="float" id="191341" name="ymin" value="0" /><Parameter format="float" id="191342" name="ymax" value="4000" /><Parameter format="int" id="191343" name="save" value="1" /><Parameter format="str" id="191344" name="figpath" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018" /></Operation></ProcUnit><ProcUnit datatype="Spectra" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="IncohInt" priority="2" type="other"><Parameter format="float" id="191221" name="n" value="2" /></Operation><Operation id="19123" name="removeDC" priority="3" type="self" /></ProcUnit></Project> No newline at end of file
@@ -1,152 +1,153
1 1 #!/usr/bin/env python
2 2 import os, sys
3 3
4 4 # path = os.path.dirname(os.getcwd())
5 5 # path = os.path.join(path, 'source')
6 6 # sys.path.insert(0, '../')
7 7
8 8 from schainpy.controller import Project
9 9
10 10 xmin = '15.5'
11 11 xmax = '24'
12 12
13 13
14 14 desc = "ProcBLTR Test"
15 15 filename = "ProcBLTR.xml"
16 16 figpath = '/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018'
17 17
18 18 controllerObj = Project()
19 19
20 20
21 21 controllerObj.setup(id='191', name='test01', description=desc)
22 22
23 23 readUnitConfObj = controllerObj.addReadUnit(datatype='BLTRSpectraReader',
24 24 #path='/media/erick/6F60F7113095A154/BLTR/',
25 25 #path='/data/BLTR',
26 path='/home/erick/Documents/Data/BLTR_Data/fdt/',
26 path = '/data/BLTR/2018/enero',
27 #path='/home/erick/Documents/Data/BLTR_Data/fdt/',
27 28 endDate='2016/11/01',
28 29 startTime='0:00:00',
29 30 startDate='2016/11/01',
30 31 endTime='23:59:59',
31 32
32 33
33 34 online=0,
34 35 walk=0,
35 36 ReadMode='1')
36 37 # expLabel='')
37 38
38 39 # opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
39 40
40 41 procUnitConfObj1 = controllerObj.addProcUnit(datatype='Spectra', inputId=readUnitConfObj.getId())
41 42
42 43
43 44 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
44 45 opObj11.addParameter(name='n', value='2', format='float')
45 46
46 47 opObj10 = procUnitConfObj1.addOperation(name='removeDC')
47 48 #opObj10 = procUnitConfObj1.addOperation(name='removeInterference2')
48 49
49 50 # opObj10 = procUnitConfObj1.addOperation(name='calcMag')
50 51
51 52 # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
52 53 # opObj11.addParameter(name='id', value='21', format='int')
53 54 # opObj11.addParameter(name='wintitle', value='SpectraCutPlot', format='str')
54 55 # opObj11.addParameter(name='xaxis', value='frequency', format='str')
55 56 # opObj11.addParameter(name='colormap', value='winter', format='str')
56 57 # opObj11.addParameter(name='xmin', value='-0.005', format='float')
57 58 # opObj11.addParameter(name='xmax', value='0.005', format='float')
58 59 # #opObj10 = procUnitConfObj1.addOperation(name='selectChannels')
59 60 # #opObj10.addParameter(name='channelList', value='0,1', format='intlist')
60 61 # opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
61 62 # opObj11.addParameter(name='id', value='21', format='int')
62 63 # opObj11.addParameter(name='wintitle', value='SpectraPlot', format='str')
63 64 # #opObj11.addParameter(name='xaxis', value='Velocity', format='str')
64 65
65 66 # opObj11.addParameter(name='xaxis', value='velocity', format='str')
66 67 # opObj11.addParameter(name='xmin', value='-0.005', format='float')
67 68 # opObj11.addParameter(name='xmax', value='0.005', format='float')
68 69
69 70 # opObj11.addParameter(name='ymin', value='225', format='float')
70 71 # opObj11.addParameter(name='ymax', value='3000', format='float')
71 72 # opObj11.addParameter(name='zmin', value='-100', format='int')
72 73 # opObj11.addParameter(name='zmax', value='-65', format='int')
73 74
74 75 # opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
75 76 # opObj11.addParameter(name='id', value='10', format='int')
76 77 # opObj11.addParameter(name='wintitle', value='RTI', format='str')
77 78 # opObj11.addParameter(name='ymin', value='0', format='float')
78 79 # opObj11.addParameter(name='ymax', value='4000', format='float')
79 80 # #opObj11.addParameter(name='zmin', value='-100', format='int')
80 81 # #opObj11.addParameter(name='zmax', value='-70', format='int')
81 82 # opObj11.addParameter(name='zmin', value='-90', format='int')
82 83 # opObj11.addParameter(name='zmax', value='-40', format='int')
83 84 # opObj11.addParameter(name='showprofile', value='1', format='int')
84 85 # opObj11.addParameter(name='timerange', value=str(2*60*60), format='int')
85 86
86 87 # opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other')
87 88 # procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
88 89 # opObj11.addParameter(name='id', value='2005', format='int')
89 90 # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str')
90 91 # # opObj11.addParameter(name='exp_code', value='13', format='int')
91 92 # opObj11.addParameter(name='xaxis', value='Velocity', format='str')
92 93 #opObj11.addParameter(name='xmin', value='-10', format='float')
93 94 #opObj11.addParameter(name='xmax', value='10', format='float')
94 95 #opObj11.addParameter(name='ymin', value='225', format='float')
95 96 #opObj11.addParameter(name='ymax', value='3000', format='float')
96 97 #opObj11.addParameter(name='phase_min', value='-4', format='int')
97 98 #opObj11.addParameter(name='phase_max', value='4', format='int')
98 99
99 100 # procUnitConfObj2 = controllerObj.addProcUnit(datatype='CorrelationProc', inputId=procUnitConfObj1.getId())
100 101 # procUnitConfObj2.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
101 102
102 103 procUnitConfObj2 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj1.getId())
103 104 #opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other')
104 105 opObj22 = procUnitConfObj2.addOperation(name='FullSpectralAnalysis', optype='other')
105 106 opObj22.addParameter(name='SNRlimit', value='2', format='float')
106 107
107 108 opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
108 109 opObj22.addParameter(name='id', value='4', format='int')
109 110 opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
110 111 #opObj22.addParameter(name='save', value='1', format='bool')
111 112 # opObj22.addParameter(name='figpath', value = '/home/erick/Pictures', format='str')
112 113
113 114 opObj22.addParameter(name='zmin', value='-20', format='int')
114 115 opObj22.addParameter(name='zmax', value='20', format='int')
115 116 opObj22.addParameter(name='zmin_ver', value='-300', format='float')
116 117 opObj22.addParameter(name='zmax_ver', value='300', format='float')
117 118 opObj22.addParameter(name='SNRmin', value='-5', format='int')
118 119 opObj22.addParameter(name='SNRmax', value='30', format='int')
119 120 # opObj22.addParameter(name='SNRthresh', value='-3.5', format='float')
120 opObj22.addParameter(name='xmin', value='0', format='float')
121 opObj22.addParameter(name='xmax', value='24', format='float')
122 opObj22.addParameter(name='ymin', value='225', format='float')
123 #opObj22.addParameter(name='ymax', value='2000', format='float')
121 opObj22.addParameter(name='xmin', value='4', format='float')
122 opObj22.addParameter(name='xmax', value='9', format='float')
123 opObj22.addParameter(name='ymin', value='0', format='float')
124 opObj22.addParameter(name='ymax', value='4000', format='float')
124 125 opObj22.addParameter(name='save', value='1', format='int')
125 126 opObj22.addParameter(name='figpath', value=figpath, format='str')
126 127
127 128
128 129 # opObj11.addParameter(name='pairlist', value='(1,0),(0,2),(1,2)', format='pairsList')
129 130 #opObj10 = procUnitConfObj1.addOperation(name='selectHeights')
130 131 #opObj10.addParameter(name='minHei', value='225', format='float')
131 132 #opObj10.addParameter(name='maxHei', value='1000', format='float')
132 133
133 134 # opObj11 = procUnitConfObj1.addOperation(name='CoherenceMap', optype='other')
134 135 # opObj11.addParameter(name='id', value='102', format='int')
135 136 # opObj11.addParameter(name='wintitle', value='Coherence', format='str')
136 137 # opObj11.addParameter(name='ymin', value='225', format='float')
137 138 # opObj11.addParameter(name='ymax', value='4000', format='float')
138 139
139 140 # opObj11.addParameter(name='phase_cmap', value='jet', format='str')
140 141 # opObj11.addParameter(name='xmin', value='8.5', format='float')
141 142 # opObj11.addParameter(name='xmax', value='9.5', format='float')
142 143 # opObj11.addParameter(name='figpath', value=figpath, format='str')
143 144 # opObj11.addParameter(name='save', value=1, format='bool')
144 145 # opObj11.addParameter(name='pairsList', value='(1,0),(3,2)', format='pairsList')
145 146
146 147 # opObj12 = procUnitConfObj1.addOperation(name='PublishData', optype='other')
147 148 # opObj12.addParameter(name='zeromq', value=1, format='int')
148 149 # opObj12.addParameter(name='verbose', value=0, format='bool')
149 150 # opObj12.addParameter(name='server', value='erick2', format='str')
150 151 controllerObj.start()
151 152
152 153
General Comments 0
You need to be logged in to leave comments. Login now