##// END OF EJS Templates
--
joabAM -
r1712:562a60b2dcd6
parent child
Show More
@@ -1,780 +1,779
1 1 ''''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @upgrade: 2021, Joab Apaza
7 7
8 8 '''
9 9
10 10 import os
11 11 import sys
12 12 import glob
13 13 import fnmatch
14 14 import datetime
15 15 import time
16 16 import re
17 17 import h5py
18 18 import numpy
19 19
20 20 try:
21 21 from gevent import sleep
22 22 except:
23 23 from time import sleep
24 24
25 25 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader,ProcessingHeader
26 26 from schainpy.model.data.jrodata import Voltage
27 27 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
28 28 from numpy import imag
29 29 from schainpy.utils import log
30 30
31 31
32 32 class AMISRReader(ProcessingUnit):
33 33 '''
34 34 classdocs
35 35 '''
36 36
37 37 def __init__(self):
38 38 '''
39 39 Constructor
40 40 '''
41 41
42 42 ProcessingUnit.__init__(self)
43 43
44 44 self.set = None
45 45 self.subset = None
46 46 self.extension_file = '.h5'
47 47 self.dtc_str = 'dtc'
48 48 self.dtc_id = 0
49 49 self.status = True
50 50 self.isConfig = False
51 51 self.dirnameList = []
52 52 self.filenameList = []
53 53 self.fileIndex = None
54 54 self.flagNoMoreFiles = False
55 55 self.flagIsNewFile = 0
56 56 self.filename = ''
57 57 self.amisrFilePointer = None
58 58
59 59 self.beamCodeMap = None
60 60 self.azimuthList = []
61 61 self.elevationList = []
62 62 self.dataShape = None
63 63 self.flag_old_beams = False
64 64
65 65 self.flagAsync = False #Use when the experiment has no syncronization
66 66 self.shiftChannels = 0
67 67 self.profileIndex = 0
68 68
69 69
70 70 self.beamCodeByFrame = None
71 71 self.radacTimeByFrame = None
72 72
73 73 self.dataset = None
74 74
75 75 self.__firstFile = True
76 76
77 77 self.buffer = None
78 78
79 79 self.timezone = 'ut'
80 80
81 81 self.__waitForNewFile = 20
82 82 self.__filename_online = None
83 83 #Is really necessary create the output object in the initializer
84 84 self.dataOut = Voltage()
85 85 self.dataOut.error=False
86 86 self.margin_days = 1
87 87 self.flag_ignoreFiles = False #to activate the ignoring Files flag
88 88 self.flag_standby = False # just keep waiting, use when ignoring files
89 89 self.ignStartDateTime=None
90 90 self.ignEndDateTime=None
91 91
92 92 def setup(self,path=None,
93 93 startDate=None,
94 94 endDate=None,
95 95 startTime=None,
96 96 endTime=None,
97 97 walk=True,
98 98 timezone='ut',
99 99 all=0,
100 100 code = 1,
101 101 nCode = 1,
102 102 nBaud = 0,
103 103 nOsamp = 0,
104 104 online=False,
105 105 old_beams=False,
106 106 margin_days=1,
107 107 nFFT = None,
108 108 nChannels = None,
109 109 ignStartDate=None,
110 110 ignEndDate=None,
111 111 ignStartTime=None,
112 112 ignEndTime=None,
113 113 syncronization=True,
114 114 shiftChannels=0
115 115 ):
116 116
117 117
118 118
119 119 self.timezone = timezone
120 120 self.all = all
121 121 self.online = online
122 122 self.flag_old_beams = old_beams
123 123 self.code = code
124 124 self.nCode = int(nCode)
125 125 self.nBaud = int(nBaud)
126 126 self.nOsamp = int(nOsamp)
127 127 self.margin_days = margin_days
128 128 self.__sampleRate = None
129 129 self.flagAsync = not syncronization
130 130 self.shiftChannels = shiftChannels
131 131 self.nFFT = nFFT
132 132 self.nChannels = nChannels
133 133 if ignStartTime!=None and ignEndTime!=None:
134 134 if ignStartDate!=None and ignEndDate!=None:
135 135 self.ignStartDateTime=datetime.datetime.combine(ignStartDate,ignStartTime)
136 136 self.ignEndDateTime=datetime.datetime.combine(ignEndDate,ignEndTime)
137 137 else:
138 138 self.ignStartDateTime=datetime.datetime.combine(startDate,ignStartTime)
139 139 self.ignEndDateTime=datetime.datetime.combine(endDate,ignEndTime)
140 140 self.flag_ignoreFiles = True
141 141
142 142 #self.findFiles()
143 143 if not(online):
144 144 #Busqueda de archivos offline
145 145 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk,)
146 146 else:
147 147 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
148 148
149 149 if not(self.filenameList):
150 150 raise schainpy.admin.SchainWarning("There is no files into the folder: %s"%(path))
151 151 #sys.exit(0)
152 152 self.dataOut.error = True
153 153
154 154 self.fileIndex = 0
155 155
156 156 self.readNextFile(online)
157 157
158 158 '''
159 159 Add code
160 160 '''
161 161 self.isConfig = True
162 162 # print("Setup Done")
163 163 pass
164 164
165 165
166 166 def readAMISRHeader(self,fp):
167 167
168 168 if self.isConfig and (not self.flagNoMoreFiles):
169 169 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
170 170 if self.dataShape != newShape and newShape != None and not self.flag_standby:
171 171 raise schainpy.admin.SchainError("NEW FILE HAS A DIFFERENT SHAPE: ")
172 172 print(self.dataShape,newShape,"\n")
173 173 return 0
174 174 else:
175 175 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
176 176
177 177
178 178 header = 'Raw11/Data/RadacHeader'
179 179 if self.nChannels == None:
180 180 expFile = fp['Setup/Experimentfile'][()].decode()
181 181 linesExp = expFile.split("\n")
182 182 a = [line for line in linesExp if "nbeamcodes" in line]
183 183 self.nChannels = int(a[0][11:])
184 184
185 185 if not self.flagAsync: #for experiments with no syncronization
186 186 self.shiftChannels = 0
187 187
188 188
189 189
190 190 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
191 191
192 192
193 193 if (self.startDate > datetime.date(2021, 7, 15)) or self.flag_old_beams: #Se cambiΓ³ la forma de extracciΓ³n de Apuntes el 17 o forzar con flag de reorganizaciΓ³n
194 194 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
195 195 self.trueBeams = self.beamcodeFile.split("\n")
196 196 self.trueBeams.pop()#remove last
197 197 if self.nFFT == None:
198 198 log.error("FFT or number of repetitions per channels is needed",self.name)
199 199 beams_idx = [k*self.nFFT for k in range(self.nChannels)]
200 200 beams = [self.trueBeams[b] for b in beams_idx]
201 201 self.beamCode = [int(x, 16) for x in beams]
202 202
203 203 if(self.flagAsync and self.shiftChannels == 0):
204 204 initBeam = self.beamCodeByPulse[0, 0]
205 205 self.shiftChannels = numpy.argwhere(self.beamCode ==initBeam)[0,0]
206 206
207 207 else:
208 208 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
209 209 self.beamCode = _beamCode[0,:]
210 210
211 211
212 212
213 213
214 214 if self.beamCodeMap == None:
215 215 self.beamCodeMap = fp['Setup/BeamcodeMap']
216 216 for beam in self.beamCode:
217 217 beamAziElev = numpy.where(self.beamCodeMap[:,0]==beam)
218 218 beamAziElev = beamAziElev[0].squeeze()
219 219 self.azimuthList.append(self.beamCodeMap[beamAziElev,1])
220 220 self.elevationList.append(self.beamCodeMap[beamAziElev,2])
221 221 #print("Beamssss: ",self.beamCodeMap[beamAziElev,1],self.beamCodeMap[beamAziElev,2])
222 222 #print(self.beamCode)
223 223 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
224 224 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
225 225 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
226 226 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
227 227 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
228 228 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
229 229 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
230 230 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
231 231 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
232 232 self.frequency = fp.get('Rx/Frequency')
233 233 txAus = fp.get('Raw11/Data/Pulsewidth') #seconds
234 234 self.baud = fp.get('Raw11/Data/TxBaud')
235 235 sampleRate = fp.get('Rx/SampleRate')
236 236 self.__sampleRate = sampleRate[()]
237 237 self.nblocks = self.pulseCount.shape[0] #nblocks
238 238 self.profPerBlockRAW = self.pulseCount.shape[1] #profiles per block in raw data
239 239 self.nprofiles = self.pulseCount.shape[1] #nprofile
240 240 #self.nsa = self.nsamplesPulse[0,0] #ngates
241 241 self.nsa = len(self.rangeFromFile[0])
242 242 self.nchannels = len(self.beamCode)
243 243 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
244 244 #print("IPPS secs: ",self.ippSeconds)
245 245 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
246 246 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
247 247
248 248 #filling radar controller header parameters
249 249 self.__ippKm = self.ippSeconds *.15*1e6 # in km
250 250 #self.__txA = txAus[()]*.15 #(ipp[us]*.15km/1us) in km
251 251 self.__txA = txAus[()] #seconds
252 252 self.__txAKm = self.__txA*1e6*.15
253 253 self.__txB = 0
254 254 nWindows=1
255 255 self.__nSamples = self.nsa
256 256 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
257 257 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
258 258 #print("amisr-ipp:",self.ippSeconds, self.__ippKm)
259 259 #for now until understand why the code saved is different (code included even though code not in tuf file)
260 260 #self.__codeType = 0
261 261 # self.__nCode = None
262 262 # self.__nBaud = None
263 263 self.__code = self.code
264 264 self.__codeType = 0
265 265 if self.code != None:
266 266 self.__codeType = 1
267 267 self.__nCode = self.nCode
268 268 self.__nBaud = self.nBaud
269 269 self.__baudTX = self.__txA/(self.nBaud)
270 270 #self.__code = 0
271 271
272 272 #filling system header parameters
273 273 self.__nSamples = self.nsa
274 274 self.newProfiles = self.nprofiles/self.nchannels
275 275 self.__channelList = [n for n in range(self.nchannels)]
276 276
277 277 self.__frequency = self.frequency[0][0]
278 278
279 279
280 280 return 1
281 281
282 282
283 283 def createBuffers(self):
284 284
285 285 pass
286 286
287 287 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
288 288 self.path = path
289 289 self.startDate = startDate
290 290 self.endDate = endDate
291 291 self.startTime = startTime
292 292 self.endTime = endTime
293 293 self.walk = walk
294 294
295 295
296 296 def __checkPath(self):
297 297 if os.path.exists(self.path):
298 298 self.status = 1
299 299 else:
300 300 self.status = 0
301 301 print('Path:%s does not exists'%self.path)
302 302
303 303 return
304 304
305 305
306 306 def __selDates(self, amisr_dirname_format):
307 307 try:
308 308 year = int(amisr_dirname_format[0:4])
309 309 month = int(amisr_dirname_format[4:6])
310 310 dom = int(amisr_dirname_format[6:8])
311 311 thisDate = datetime.date(year,month,dom)
312 312 #margen de un dΓ­a extra, igual luego se filtra for fecha y hora
313 313 if (thisDate>=(self.startDate - datetime.timedelta(days=self.margin_days)) and thisDate <= (self.endDate)+ datetime.timedelta(days=1)):
314 314 return amisr_dirname_format
315 315 except:
316 316 return None
317 317
318 318
319 319 def __findDataForDates(self,online=False):
320 320
321 321 if not(self.status):
322 322 return None
323 323
324 324 pat = '\d+.\d+'
325 325 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
326 326 dirnameList = [x for x in dirnameList if x!=None]
327 327 dirnameList = [x.string for x in dirnameList]
328 328 if not(online):
329 329 dirnameList = [self.__selDates(x) for x in dirnameList]
330 330 dirnameList = [x for x in dirnameList if x!=None]
331 331 if len(dirnameList)>0:
332 332 self.status = 1
333 333 self.dirnameList = dirnameList
334 334 self.dirnameList.sort()
335 335 else:
336 336 self.status = 0
337 337 return None
338 338
339 339 def __getTimeFromData(self):
340 340 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
341 341 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
342 342
343 343 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
344 344 print('........................................')
345 345 filter_filenameList = []
346 346 self.filenameList.sort()
347 347 total_files = len(self.filenameList)
348 348 #for i in range(len(self.filenameList)-1):
349 349 for i in range(total_files):
350 350 filename = self.filenameList[i]
351 351 #print("file-> ",filename)
352 352 try:
353 353 fp = h5py.File(filename,'r')
354 354 time_str = fp.get('Time/RadacTimeString')
355 355
356 356 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
357 357 #startDateTimeStr_File = "2019-12-16 09:21:11"
358 358 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
359 359 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
360 360
361 361 #endDateTimeStr_File = "2019-12-16 11:10:11"
362 362 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
363 363 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
364 364 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
365 365
366 366 fp.close()
367 367
368 368 #print("check time", startDateTime_File)
369 369 if self.timezone == 'lt':
370 370 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
371 371 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
372 372 if (startDateTime_File >=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
373 373 filter_filenameList.append(filename)
374 374
375 375 if (startDateTime_File>endDateTime_Reader):
376 376 break
377 377 except Exception as e:
378 378 log.warning("Error opening file {} -> {}".format(os.path.split(filename)[1],e))
379 379
380 380 filter_filenameList.sort()
381 381 self.filenameList = filter_filenameList
382 382
383 383 return 1
384 384
385 385 def __filterByGlob1(self, dirName):
386 386 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
387 387 filter_files.sort()
388 388 filterDict = {}
389 389 filterDict.setdefault(dirName)
390 390 filterDict[dirName] = filter_files
391 391 return filterDict
392 392
393 393 def __getFilenameList(self, fileListInKeys, dirList):
394 394 for value in fileListInKeys:
395 395 dirName = list(value.keys())[0]
396 396 for file in value[dirName]:
397 397 filename = os.path.join(dirName, file)
398 398 self.filenameList.append(filename)
399 399
400 400
401 401 def __selectDataForTimes(self, online=False):
402 402 #aun no esta implementado el filtro for tiempo-> implementado en readNextFile
403 403 if not(self.status):
404 404 return None
405 405
406 406 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
407 407 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
408 408 self.__getFilenameList(fileListInKeys, dirList)
409 409 if not(online):
410 410 #filtro por tiempo
411 411 if not(self.all):
412 412 self.__getTimeFromData()
413 413
414 414 if len(self.filenameList)>0:
415 415 self.status = 1
416 416 self.filenameList.sort()
417 417 else:
418 418 self.status = 0
419 419 return None
420 420
421 421 else:
422 422 #get the last file - 1
423 423 self.filenameList = [self.filenameList[-2]]
424 424 new_dirnameList = []
425 425 for dirname in self.dirnameList:
426 426 junk = numpy.array([dirname in x for x in self.filenameList])
427 427 junk_sum = junk.sum()
428 428 if junk_sum > 0:
429 429 new_dirnameList.append(dirname)
430 430 self.dirnameList = new_dirnameList
431 431 return 1
432 432
433 433 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
434 434 endTime=datetime.time(23,59,59),walk=True):
435 435
436 436 if endDate ==None:
437 437 startDate = datetime.datetime.utcnow().date()
438 438 endDate = datetime.datetime.utcnow().date()
439 439
440 440 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
441 441
442 442 self.__checkPath()
443 443
444 444 self.__findDataForDates(online=True)
445 445
446 446 self.dirnameList = [self.dirnameList[-1]]
447 447
448 448 self.__selectDataForTimes(online=True)
449 449
450 450 return
451 451
452 452
453 453 def searchFilesOffLine(self,
454 454 path,
455 455 startDate,
456 456 endDate,
457 457 startTime=datetime.time(0,0,0),
458 458 endTime=datetime.time(23,59,59),
459 459 walk=True):
460 460
461 461 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
462 462
463 463 self.__checkPath()
464 464
465 465 self.__findDataForDates()
466 466
467 467 self.__selectDataForTimes()
468 468
469 469 for i in range(len(self.filenameList)):
470 470 print("%s" %(self.filenameList[i]))
471 471
472 472 return
473 473
474 474 def __setNextFileOffline(self):
475 475
476 476 try:
477 477 self.filename = self.filenameList[self.fileIndex]
478 478 self.amisrFilePointer = h5py.File(self.filename,'r')
479 479 self.fileIndex += 1
480 480 except:
481 481 self.flagNoMoreFiles = 1
482 482 raise schainpy.admin.SchainError('No more files to read')
483 483 return 0
484 484
485 485 self.flagIsNewFile = 1
486 486 print("Setting the file: %s"%self.filename)
487 487
488 488 return 1
489 489
490 490
491 491 def __setNextFileOnline(self):
492 492 filename = self.filenameList[0]
493 493 if self.__filename_online != None:
494 494 self.__selectDataForTimes(online=True)
495 495 filename = self.filenameList[0]
496 496 wait = 0
497 497 self.__waitForNewFile=300 ## DEBUG:
498 498 while self.__filename_online == filename:
499 499 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
500 500 if wait == 5:
501 501 self.flagNoMoreFiles = 1
502 502 return 0
503 503 sleep(self.__waitForNewFile)
504 504 self.__selectDataForTimes(online=True)
505 505 filename = self.filenameList[0]
506 506 wait += 1
507 507
508 508 self.__filename_online = filename
509 509
510 510 self.amisrFilePointer = h5py.File(filename,'r')
511 511 self.flagIsNewFile = 1
512 512 self.filename = filename
513 513 print("Setting the file: %s"%self.filename)
514 514 return 1
515 515
516 516
517 517 def readData(self):
518 518 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
519 519 re = buffer[:,:,:,0]
520 520 im = buffer[:,:,:,1]
521 521 dataset = re + im*1j
522 522
523 523 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
524 524 timeset = self.radacTime[:,0]
525 525
526 526 return dataset,timeset
527 527
528 528 def reshapeData(self):
529 529 #print(self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa)
530 530 channels = self.beamCodeByPulse[0,:]
531 531 nchan = self.nchannels
532 532 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
533 533 nblocks = self.nblocks
534 534 nsamples = self.nsa
535 535 #print("Channels: ",self.nChannels)
536 print("dataset: ", self.dataset.shape)
537 536 #Dimensions : nChannels, nProfiles, nSamples
538 537 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
539 538 ############################################
540 539 profPerCH = int(self.profPerBlockRAW / (self.nFFT* self.nChannels))
541 540 #profPerCH = int(self.profPerBlockRAW / self.nChannels)
542 541 for thisChannel in range(nchan):
543 542
544 543 ich = thisChannel
545 544
546 545 idx_ch = [self.nFFT*(ich + nchan*k) for k in range(profPerCH)]
547 546 #print(idx_ch)
548 547 if self.nFFT > 1:
549 548 aux = [numpy.arange(i, i+self.nFFT) for i in idx_ch]
550 549 idx_ch = None
551 550 idx_ch =aux
552 551 idx_ch = numpy.array(idx_ch, dtype=int).flatten()
553 552 else:
554 553 idx_ch = numpy.array(idx_ch, dtype=int)
555 554
556 555 #print(ich,profPerCH,idx_ch)
557 556 #print(numpy.where(channels==self.beamCode[ich])[0])
558 557 #new_block[:,ich,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[ich])[0],:]
559 558 new_block[:,ich,:,:] = self.dataset[:,idx_ch,:]
560 559
561 560 new_block = numpy.transpose(new_block, (1,0,2,3))
562 561 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
563 562 if self.flagAsync:
564 563 new_block = numpy.roll(new_block, self.shiftChannels, axis=0)
565 564 return new_block
566 565
567 566 def updateIndexes(self):
568 567
569 568 pass
570 569
571 570 def fillJROHeader(self):
572 571
573 572 #fill radar controller header
574 573
575 574 #fill system header
576 575 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
577 576 nProfiles=self.newProfiles,
578 577 nChannels=len(self.__channelList),
579 578 adcResolution=14,
580 579 pciDioBusWidth=32)
581 580
582 581 self.dataOut.type = "Voltage"
583 582 self.dataOut.data = None
584 583 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
585 584 # self.dataOut.nChannels = 0
586 585
587 586 # self.dataOut.nHeights = 0
588 587
589 588 self.dataOut.nProfiles = self.newProfiles*self.nblocks
590 589 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
591 590 ranges = numpy.reshape(self.rangeFromFile[()],(-1))
592 591 self.dataOut.heightList = ranges/1000.0 #km
593 592 self.dataOut.channelList = self.__channelList
594 593
595 594 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
596 595
597 596 # self.dataOut.channelIndexList = None
598 597
599 598
600 599 # #self.dataOut.azimuthList = numpy.roll( numpy.array(self.azimuthList) ,self.shiftChannels)
601 600 # #self.dataOut.elevationList = numpy.roll(numpy.array(self.elevationList) ,self.shiftChannels)
602 601 # #self.dataOut.codeList = numpy.roll(numpy.array(self.beamCode), self.shiftChannels)
603 602
604 603 self.dataOut.azimuthList = self.azimuthList
605 604 self.dataOut.elevationList = self.elevationList
606 605 self.dataOut.codeList = self.beamCode
607 606
608 607
609 608
610 609 #print(self.dataOut.elevationList)
611 610 self.dataOut.flagNoData = True
612 611
613 612 #Set to TRUE if the data is discontinuous
614 613 self.dataOut.flagDiscontinuousBlock = False
615 614
616 615 self.dataOut.utctime = None
617 616
618 617 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
619 618 if self.timezone == 'lt':
620 619 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
621 620 else:
622 621 self.dataOut.timeZone = 0 #by default time is UTC
623 622
624 623 self.dataOut.dstFlag = 0
625 624 self.dataOut.errorCount = 0
626 625 self.dataOut.nCohInt = 1
627 626 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
628 627 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
629 628 self.dataOut.flagShiftFFT = False
630 629 self.dataOut.ippSeconds = self.ippSeconds
631 630 self.dataOut.ipp = self.__ippKm
632 631 self.dataOut.nCode = self.__nCode
633 632 self.dataOut.code = self.__code
634 633 self.dataOut.nBaud = self.__nBaud
635 634
636 635
637 636 self.dataOut.frequency = self.__frequency
638 637 self.dataOut.realtime = self.online
639 638
640 639 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
641 640 txA=self.__txAKm,
642 641 txB=0,
643 642 nWindows=1,
644 643 nHeights=self.__nSamples,
645 644 firstHeight=self.__firstHeight,
646 645 codeType=self.__codeType,
647 646 nCode=self.__nCode, nBaud=self.__nBaud,
648 647 code = self.__code,
649 648 nOsamp=self.nOsamp,
650 649 frequency = self.__frequency,
651 650 sampleRate= self.__sampleRate,
652 651 fClock=self.__sampleRate)
653 652
654 653
655 654 self.dataOut.radarControllerHeaderObj.heightList = ranges/1000.0 #km
656 655 self.dataOut.radarControllerHeaderObj.heightResolution = self.__deltaHeight
657 656 self.dataOut.radarControllerHeaderObj.rangeIpp = self.__ippKm #km
658 657 self.dataOut.radarControllerHeaderObj.rangeTxA = self.__txA*1e6*.15 #km
659 658 self.dataOut.radarControllerHeaderObj.nChannels = self.nchannels
660 659 self.dataOut.radarControllerHeaderObj.channelList = self.__channelList
661 660 self.dataOut.radarControllerHeaderObj.azimuthList = self.azimuthList
662 661 self.dataOut.radarControllerHeaderObj.elevationList = self.elevationList
663 662 self.dataOut.radarControllerHeaderObj.dtype = "Voltage"
664 663 self.dataOut.ippSeconds = self.ippSeconds
665 664 self.dataOut.ippFactor = self.nFFT
666 665 pass
667 666
668 667 def readNextFile(self,online=False):
669 668
670 669 if not(online):
671 670 newFile = self.__setNextFileOffline()
672 671 else:
673 672 newFile = self.__setNextFileOnline()
674 673
675 674 if not(newFile):
676 675 self.dataOut.error = True
677 676 return 0
678 677
679 678 if not self.readAMISRHeader(self.amisrFilePointer):
680 679 self.dataOut.error = True
681 680 return 0
682 681
683 682 #self.createBuffers()
684 683 self.fillJROHeader()
685 684
686 685 #self.__firstFile = False
687 686
688 687 self.dataset,self.timeset = self.readData()
689 688
690 689 if self.endDate!=None:
691 690 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
692 691 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
693 692 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
694 693 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
695 694 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
696 695 if self.timezone == 'lt':
697 696 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
698 697 if (startDateTime_File>endDateTime_Reader):
699 698 self.flag_standby = False
700 699 return 0
701 700 if self.flag_ignoreFiles and (startDateTime_File >= self.ignStartDateTime and startDateTime_File <= self.ignEndDateTime):
702 701 print("Ignoring...")
703 702 self.flag_standby = True
704 703 return 1
705 704 self.flag_standby = False
706 705
707 706 self.jrodataset = self.reshapeData()
708 707 #----self.updateIndexes()
709 708 self.profileIndex = 0
710 709
711 710 return 1
712 711
713 712
714 713 def __hasNotDataInBuffer(self):
715 714 if self.profileIndex >= (self.newProfiles*self.nblocks):
716 715 return 1
717 716 return 0
718 717
719 718
720 719 def getData(self):
721 720
722 721 if self.flagNoMoreFiles:
723 722 self.dataOut.flagNoData = True
724 723 return 0
725 724
726 725 if self.profileIndex >= (self.newProfiles*self.nblocks): #
727 726 #if self.__hasNotDataInBuffer():
728 727 if not (self.readNextFile(self.online)):
729 728 print("Profile Index break...")
730 729 return 0
731 730
732 731 if self.flag_standby: #Standby mode, if files are being ignoring, just return with no error flag
733 732 return 0
734 733
735 734 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
736 735 self.dataOut.flagNoData = True
737 736 print("No more data break...")
738 737 return 0
739 738
740 739 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
741 740
742 741 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
743 742
744 743 #print("R_t",self.timeset)
745 744
746 745 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
747 746 #verificar basic header de jro data y ver si es compatible con este valor
748 747 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
749 748 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
750 749 indexblock = self.profileIndex/self.newProfiles
751 750 #print (indexblock, indexprof)
752 751 diffUTC = 0
753 752 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
754 753
755 754 #print("utc :",indexblock," __ ",t_comp)
756 755 #print(numpy.shape(self.timeset))
757 756 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
758 757 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
759 758
760 759 self.dataOut.profileIndex = self.profileIndex
761 760 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
762 761 self.dataOut.flagNoData = False
763 762 # if indexprof == 0:
764 763 # print("kamisr: ",self.dataOut.utctime)
765 764
766 765 self.profileIndex += 1
767 766
768 767 return self.dataOut.data #retorno necesario??
769 768
770 769
771 770 def run(self, **kwargs):
772 771 '''
773 772 This method will be called many times so here you should put all your code
774 773 '''
775 774 #print("running kamisr")
776 775 if not self.isConfig:
777 776 self.setup(**kwargs)
778 777 self.isConfig = True
779 778
780 779 self.getData()
@@ -1,2306 +1,2305
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.model.data import _noise
21 21
22 22 from schainpy.utils import log
23 23 import matplotlib.pyplot as plt
24 24 #from scipy.optimize import curve_fit
25 25 from schainpy.model.io.utilsIO import getHei_index
26 26 import datetime
27 27
28 28 class SpectraProc(ProcessingUnit):
29 29
30 30 def __init__(self):
31 31
32 32 ProcessingUnit.__init__(self)
33 33
34 34 self.buffer = None
35 35 self.firstdatatime = None
36 36 self.profIndex = 0
37 37 self.dataOut = Spectra()
38 38 self.id_min = None
39 39 self.id_max = None
40 40 self.setupReq = False #Agregar a todas las unidades de proc
41 41 self.nsamplesFFT = 0
42 42
43 43 def __updateSpecFromVoltage(self):
44 44
45 45
46 46
47 47 self.dataOut.timeZone = self.dataIn.timeZone
48 48 self.dataOut.dstFlag = self.dataIn.dstFlag
49 49 self.dataOut.errorCount = self.dataIn.errorCount
50 50 self.dataOut.useLocalTime = self.dataIn.useLocalTime
51 51
52 52 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
53 53 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
54 54 self.dataOut.ippSeconds = self.dataIn.ippSeconds
55 55 self.dataOut.ipp = self.dataIn.ipp
56 56 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
57 57 self.dataOut.channelList = self.dataIn.channelList
58 58 self.dataOut.heightList = self.dataIn.heightList
59 59 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
60 60 self.dataOut.nProfiles = self.dataOut.nFFTPoints
61 61 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
62 62 self.dataOut.utctime = self.firstdatatime
63 63 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
64 64 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
65 65 self.dataOut.flagShiftFFT = False
66 66 self.dataOut.nCohInt = self.dataIn.nCohInt
67 67 self.dataOut.nIncohInt = 1
68 68 self.dataOut.deltaHeight = self.dataIn.deltaHeight
69 69 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
70 70 self.dataOut.frequency = self.dataIn.frequency
71 71 self.dataOut.realtime = self.dataIn.realtime
72 72 self.dataOut.azimuth = self.dataIn.azimuth
73 73 self.dataOut.zenith = self.dataIn.zenith
74 74 self.dataOut.codeList = self.dataIn.codeList
75 75 self.dataOut.azimuthList = self.dataIn.azimuthList
76 76 self.dataOut.elevationList = self.dataIn.elevationList
77 77 self.dataOut.code = self.dataIn.code
78 78 self.dataOut.nCode = self.dataIn.nCode
79 79 self.dataOut.flagProfilesByRange = self.dataIn.flagProfilesByRange
80 80 self.dataOut.nProfilesByRange = self.dataIn.nProfilesByRange
81 81
82 82
83 83 def __getFft(self):
84 84 # print("fft donw")
85 85 """
86 86 Convierte valores de Voltaje a Spectra
87 87
88 88 Affected:
89 89 self.dataOut.data_spc
90 90 self.dataOut.data_cspc
91 91 self.dataOut.data_dc
92 92 self.dataOut.heightList
93 93 self.profIndex
94 94 self.buffer
95 95 self.dataOut.flagNoData
96 96 """
97 97 fft_volt = numpy.fft.fft(
98 98 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
99 99 fft_volt = fft_volt.astype(numpy.dtype('complex'))
100 100 dc = fft_volt[:, 0, :]
101 101
102 102 # calculo de self-spectra
103 103 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
104 104 spc = fft_volt * numpy.conjugate(fft_volt)
105 105 spc = spc.real
106 106
107 107 blocksize = 0
108 108 blocksize += dc.size
109 109 blocksize += spc.size
110 110
111 111 cspc = None
112 112 pairIndex = 0
113 113 if self.dataOut.pairsList != None:
114 114 # calculo de cross-spectra
115 115 cspc = numpy.zeros(
116 116 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
117 117 for pair in self.dataOut.pairsList:
118 118 if pair[0] not in self.dataOut.channelList:
119 119 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
120 120 str(pair), str(self.dataOut.channelList)))
121 121 if pair[1] not in self.dataOut.channelList:
122 122 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
123 123 str(pair), str(self.dataOut.channelList)))
124 124
125 125 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
126 126 numpy.conjugate(fft_volt[pair[1], :, :])
127 127 pairIndex += 1
128 128 blocksize += cspc.size
129 129
130 130 self.dataOut.data_spc = spc
131 131 self.dataOut.data_cspc = cspc
132 132 self.dataOut.data_dc = dc
133 133 self.dataOut.blockSize = blocksize
134 134 self.dataOut.flagShiftFFT = False
135 135
136 136 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, zeroPad=False, zeroPoints=0):
137 137
138 138
139 139 try:
140 140 type = self.dataIn.type.decode("utf-8")
141 141 self.dataIn.type = type
142 142 except:
143 143 pass
144 144 if self.dataIn.type == "Spectra":
145 145
146 146 try:
147 147 self.dataOut.copy(self.dataIn)
148 148 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
149 149 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
150 150 self.dataOut.nProfiles = self.dataOut.nFFTPoints
151 151 #self.dataOut.nHeights = len(self.dataOut.heightList)
152 152 except Exception as e:
153 153 print("Error dataIn ",e)
154 154
155 155
156 156
157 157 if shift_fft:
158 158 #desplaza a la derecha en el eje 2 determinadas posiciones
159 159 shift = int(self.dataOut.nFFTPoints/2)
160 160 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
161 161
162 162 if self.dataOut.data_cspc is not None:
163 163 #desplaza a la derecha en el eje 2 determinadas posiciones
164 164 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
165 165 if pairsList:
166 166 self.__selectPairs(pairsList)
167 167
168 168
169 169 elif self.dataIn.type == "Voltage":
170 170
171 171 self.dataOut.flagNoData = True
172 172 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
173 173 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
174 174 if nFFTPoints == None:
175 175 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
176 176
177 177 if nProfiles == None:
178 178 nProfiles = nFFTPoints
179 179
180 180 #if ippFactor == None:
181 181 # self.dataOut.ippFactor = 1
182 182 if ippFactor == None:
183 183 self.dataOut.ippFactor = self.dataIn.ippFactor
184 184 else:
185 185 self.dataOut.ippFactor = ippFactor
186 186
187 187 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
188 188 if self.buffer is None:
189 189 if not zeroPad:
190 190 self.buffer = numpy.zeros((self.dataIn.nChannels,
191 191 nProfiles,
192 192 self.dataIn.nHeights),
193 193 dtype='complex')
194 194 zeroPoints = 0
195 195 else:
196 196 self.buffer = numpy.zeros((self.dataIn.nChannels,
197 197 nFFTPoints+int(zeroPoints),
198 198 self.dataIn.nHeights),
199 199 dtype='complex')
200 200
201 201 self.dataOut.nFFTPoints = nFFTPoints + int(zeroPoints)
202 202
203 203 if self.dataIn.flagDataAsBlock:
204 204 nVoltProfiles = self.dataIn.data.shape[1]
205 205 zeroPoints = 0
206 206 if nVoltProfiles == nProfiles or zeroPad:
207 207 self.buffer = self.dataIn.data.copy()
208 208 self.profIndex = nVoltProfiles
209 209
210 210 elif nVoltProfiles < nProfiles:
211 211
212 212 if self.profIndex == 0:
213 213 self.id_min = 0
214 214 self.id_max = nVoltProfiles
215 215
216 216 self.buffer[:, self.id_min:self.id_max,
217 217 :] = self.dataIn.data
218 218 self.profIndex += nVoltProfiles
219 219 self.id_min += nVoltProfiles
220 220 self.id_max += nVoltProfiles
221 221 else:
222 222 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
223 223 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
224 224 self.dataOut.flagNoData = True
225 225 else:
226 226 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
227 227 self.profIndex += 1
228 228
229 229 if self.firstdatatime == None:
230 230 self.firstdatatime = self.dataIn.utctime
231 231
232 232 if self.profIndex == nProfiles or (zeroPad and zeroPoints==0):
233 233
234 234 self.__updateSpecFromVoltage()
235 235
236 236 if pairsList == None:
237 237 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
238 238 else:
239 239 self.dataOut.pairsList = pairsList
240 240 self.__getFft()
241 241 self.dataOut.flagNoData = False
242 242 self.firstdatatime = None
243 243 self.nsamplesFFT = self.profIndex
244 244 self.profIndex = 0
245 245
246 246 #update Processing Header:
247 247 self.dataOut.processingHeaderObj.dtype = "Spectra"
248 248 self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
249 249 self.dataOut.processingHeaderObj.nSamplesFFT = self.nsamplesFFT
250 250 self.dataOut.processingHeaderObj.nIncohInt = 1
251 251
252 252
253 253 elif self.dataIn.type == "Parameters":
254 254
255 255 self.dataOut.data_spc = self.dataIn.data_spc
256 256 self.dataOut.data_cspc = self.dataIn.data_cspc
257 257 self.dataOut.data_outlier = self.dataIn.data_outlier
258 258 self.dataOut.nProfiles = self.dataIn.nProfiles
259 259 self.dataOut.nIncohInt = self.dataIn.nIncohInt
260 260 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
261 261 self.dataOut.ippFactor = self.dataIn.ippFactor
262 262 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
263 263 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
264 264 self.dataOut.ProcessingHeader = self.dataIn.ProcessingHeader.copy()
265 265 self.dataOut.ippSeconds = self.dataIn.ippSeconds
266 266 self.dataOut.ipp = self.dataIn.ipp
267 267 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
268 268 #self.dataOut.spc_noise = self.dataIn.getNoise()
269 269 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
270 270 # self.dataOut.normFactor = self.dataIn.normFactor
271 271 if hasattr(self.dataIn, 'channelList'):
272 272 self.dataOut.channelList = self.dataIn.channelList
273 273 if hasattr(self.dataIn, 'pairsList'):
274 274 self.dataOut.pairsList = self.dataIn.pairsList
275 275 self.dataOut.groupList = self.dataIn.pairsList
276 276
277 277 self.dataOut.flagNoData = False
278 278
279 279 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
280 280 self.dataOut.ChanDist = self.dataIn.ChanDist
281 281 else: self.dataOut.ChanDist = None
282 282
283 283 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
284 284 # self.dataOut.VelRange = self.dataIn.VelRange
285 285 #else: self.dataOut.VelRange = None
286 286
287 287
288 288
289 289 else:
290 290 raise ValueError("The type of input object {} is not valid".format(
291 291 self.dataIn.type))
292 292
293 293
294 294
295 295
296 296 #print("spc proc Done", self.dataOut.data_spc.shape)
297 297 #print(self.dataOut.data_spc)
298 298 return
299 299
300 300 def __selectPairs(self, pairsList):
301 301
302 302 if not pairsList:
303 303 return
304 304
305 305 pairs = []
306 306 pairsIndex = []
307 307
308 308 for pair in pairsList:
309 309 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
310 310 continue
311 311 pairs.append(pair)
312 312 pairsIndex.append(pairs.index(pair))
313 313
314 314 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
315 315 self.dataOut.pairsList = pairs
316 316
317 317 return
318 318
319 319 def selectFFTs(self, minFFT, maxFFT ):
320 320 """
321 321 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
322 322 minFFT<= FFT <= maxFFT
323 323 """
324 324
325 325 if (minFFT > maxFFT):
326 326 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
327 327
328 328 if (minFFT < self.dataOut.getFreqRange()[0]):
329 329 minFFT = self.dataOut.getFreqRange()[0]
330 330
331 331 if (maxFFT > self.dataOut.getFreqRange()[-1]):
332 332 maxFFT = self.dataOut.getFreqRange()[-1]
333 333
334 334 minIndex = 0
335 335 maxIndex = 0
336 336 FFTs = self.dataOut.getFreqRange()
337 337
338 338 inda = numpy.where(FFTs >= minFFT)
339 339 indb = numpy.where(FFTs <= maxFFT)
340 340
341 341 try:
342 342 minIndex = inda[0][0]
343 343 except:
344 344 minIndex = 0
345 345
346 346 try:
347 347 maxIndex = indb[0][-1]
348 348 except:
349 349 maxIndex = len(FFTs)
350 350
351 351 self.selectFFTsByIndex(minIndex, maxIndex)
352 352
353 353 return 1
354 354
355 355 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
356 356 newheis = numpy.where(
357 357 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
358 358
359 359 if hei_ref != None:
360 360 newheis = numpy.where(self.dataOut.heightList > hei_ref)
361 361
362 362 minIndex = min(newheis[0])
363 363 maxIndex = max(newheis[0])
364 364 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
365 365 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
366 366
367 367 # determina indices
368 368 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
369 369 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
370 370 avg_dB = 10 * \
371 371 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
372 372 beacon_dB = numpy.sort(avg_dB)[-nheis:]
373 373 beacon_heiIndexList = []
374 374 for val in avg_dB.tolist():
375 375 if val >= beacon_dB[0]:
376 376 beacon_heiIndexList.append(avg_dB.tolist().index(val))
377 377
378 378 #data_spc = data_spc[:,:,beacon_heiIndexList]
379 379 data_cspc = None
380 380 if self.dataOut.data_cspc is not None:
381 381 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
382 382 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
383 383
384 384 data_dc = None
385 385 if self.dataOut.data_dc is not None:
386 386 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
387 387 #data_dc = data_dc[:,beacon_heiIndexList]
388 388
389 389 self.dataOut.data_spc = data_spc
390 390 self.dataOut.data_cspc = data_cspc
391 391 self.dataOut.data_dc = data_dc
392 392 self.dataOut.heightList = heightList
393 393 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
394 394
395 395 return 1
396 396
397 397 def selectFFTsByIndex(self, minIndex, maxIndex):
398 398 """
399 399
400 400 """
401 401
402 402 if (minIndex < 0) or (minIndex > maxIndex):
403 403 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
404 404
405 405 if (maxIndex >= self.dataOut.nProfiles):
406 406 maxIndex = self.dataOut.nProfiles-1
407 407
408 408 #Spectra
409 409 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
410 410
411 411 data_cspc = None
412 412 if self.dataOut.data_cspc is not None:
413 413 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
414 414
415 415 data_dc = None
416 416 if self.dataOut.data_dc is not None:
417 417 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
418 418
419 419 self.dataOut.data_spc = data_spc
420 420 self.dataOut.data_cspc = data_cspc
421 421 self.dataOut.data_dc = data_dc
422 422
423 423 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
424 424 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
425 425 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
426 426
427 427 return 1
428 428
429 429 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
430 430 # validacion de rango
431 431 if minHei == None:
432 432 minHei = self.dataOut.heightList[0]
433 433
434 434 if maxHei == None:
435 435 maxHei = self.dataOut.heightList[-1]
436 436
437 437 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
438 438 print('minHei: %.2f is out of the heights range' % (minHei))
439 439 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
440 440 minHei = self.dataOut.heightList[0]
441 441
442 442 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
443 443 print('maxHei: %.2f is out of the heights range' % (maxHei))
444 444 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
445 445 maxHei = self.dataOut.heightList[-1]
446 446
447 447 # validacion de velocidades
448 448 velrange = self.dataOut.getVelRange(1)
449 449
450 450 if minVel == None:
451 451 minVel = velrange[0]
452 452
453 453 if maxVel == None:
454 454 maxVel = velrange[-1]
455 455
456 456 if (minVel < velrange[0]) or (minVel > maxVel):
457 457 print('minVel: %.2f is out of the velocity range' % (minVel))
458 458 print('minVel is setting to %.2f' % (velrange[0]))
459 459 minVel = velrange[0]
460 460
461 461 if (maxVel > velrange[-1]) or (maxVel < minVel):
462 462 print('maxVel: %.2f is out of the velocity range' % (maxVel))
463 463 print('maxVel is setting to %.2f' % (velrange[-1]))
464 464 maxVel = velrange[-1]
465 465
466 466 # seleccion de indices para rango
467 467 minIndex = 0
468 468 maxIndex = 0
469 469 heights = self.dataOut.heightList
470 470
471 471 inda = numpy.where(heights >= minHei)
472 472 indb = numpy.where(heights <= maxHei)
473 473
474 474 try:
475 475 minIndex = inda[0][0]
476 476 except:
477 477 minIndex = 0
478 478
479 479 try:
480 480 maxIndex = indb[0][-1]
481 481 except:
482 482 maxIndex = len(heights)
483 483
484 484 if (minIndex < 0) or (minIndex > maxIndex):
485 485 raise ValueError("some value in (%d,%d) is not valid" % (
486 486 minIndex, maxIndex))
487 487
488 488 if (maxIndex >= self.dataOut.nHeights):
489 489 maxIndex = self.dataOut.nHeights - 1
490 490
491 491 # seleccion de indices para velocidades
492 492 indminvel = numpy.where(velrange >= minVel)
493 493 indmaxvel = numpy.where(velrange <= maxVel)
494 494 try:
495 495 minIndexVel = indminvel[0][0]
496 496 except:
497 497 minIndexVel = 0
498 498
499 499 try:
500 500 maxIndexVel = indmaxvel[0][-1]
501 501 except:
502 502 maxIndexVel = len(velrange)
503 503
504 504 # seleccion del espectro
505 505 data_spc = self.dataOut.data_spc[:,
506 506 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
507 507 # estimacion de ruido
508 508 noise = numpy.zeros(self.dataOut.nChannels)
509 509
510 510 for channel in range(self.dataOut.nChannels):
511 511 daux = data_spc[channel, :, :]
512 512 sortdata = numpy.sort(daux, axis=None)
513 513 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
514 514
515 515 self.dataOut.noise_estimation = noise.copy()
516 516
517 517 return 1
518 518
519 519 class removeDC(Operation):
520 520
521 521 def run(self, dataOut, mode=2):
522 522 self.dataOut = dataOut
523 523 jspectra = self.dataOut.data_spc
524 524 jcspectra = self.dataOut.data_cspc
525 525
526 526 num_chan = jspectra.shape[0]
527 527 num_hei = jspectra.shape[2]
528 528
529 529 if jcspectra is not None:
530 530 jcspectraExist = True
531 531 num_pairs = jcspectra.shape[0]
532 532 else:
533 533 jcspectraExist = False
534 534
535 535 freq_dc = int(jspectra.shape[1] / 2)
536 536 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
537 537 ind_vel = ind_vel.astype(int)
538 538
539 539 if ind_vel[0] < 0:
540 540 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
541 541
542 542 if mode == 1:
543 543 jspectra[:, freq_dc, :] = (
544 544 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
545 545
546 546 if jcspectraExist:
547 547 jcspectra[:, freq_dc, :] = (
548 548 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
549 549
550 550 if mode == 2:
551 551
552 552 vel = numpy.array([-2, -1, 1, 2])
553 553 xx = numpy.zeros([4, 4])
554 554
555 555 for fil in range(4):
556 556 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
557 557
558 558 xx_inv = numpy.linalg.inv(xx)
559 559 xx_aux = xx_inv[0, :]
560 560
561 561 for ich in range(num_chan):
562 562 yy = jspectra[ich, ind_vel, :]
563 563 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
564 564
565 565 junkid = jspectra[ich, freq_dc, :] <= 0
566 566 cjunkid = sum(junkid)
567 567
568 568 if cjunkid.any():
569 569 jspectra[ich, freq_dc, junkid.nonzero()] = (
570 570 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
571 571
572 572 if jcspectraExist:
573 573 for ip in range(num_pairs):
574 574 yy = jcspectra[ip, ind_vel, :]
575 575 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
576 576
577 577 self.dataOut.data_spc = jspectra
578 578 self.dataOut.data_cspc = jcspectra
579 579
580 580 return self.dataOut
581 581
582 582 class getNoiseB(Operation):
583 583
584 584 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
585 585 def __init__(self):
586 586
587 587 Operation.__init__(self)
588 588 self.isConfig = False
589 589
590 590 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
591 591
592 592 self.warnings = warnings
593 593 if minHei == None:
594 594 minHei = self.dataOut.heightList[0]
595 595
596 596 if maxHei == None:
597 597 maxHei = self.dataOut.heightList[-1]
598 598
599 599 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
600 600 if self.warnings:
601 601 print('minHei: %.2f is out of the heights range' % (minHei))
602 602 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
603 603 minHei = self.dataOut.heightList[0]
604 604
605 605 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
606 606 if self.warnings:
607 607 print('maxHei: %.2f is out of the heights range' % (maxHei))
608 608 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
609 609 maxHei = self.dataOut.heightList[-1]
610 610
611 611
612 612 #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
613 613 minIndexFFT = 0
614 614 maxIndexFFT = 0
615 615 # validacion de velocidades
616 616 indminPoint = None
617 617 indmaxPoint = None
618 618 if self.dataOut.type == 'Spectra':
619 619 if minVel == None and maxVel == None :
620 620
621 621 freqrange = self.dataOut.getFreqRange(1)
622 622
623 623 if minFreq == None:
624 624 minFreq = freqrange[0]
625 625
626 626 if maxFreq == None:
627 627 maxFreq = freqrange[-1]
628 628
629 629 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
630 630 if self.warnings:
631 631 print('minFreq: %.2f is out of the frequency range' % (minFreq))
632 632 print('minFreq is setting to %.2f' % (freqrange[0]))
633 633 minFreq = freqrange[0]
634 634
635 635 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
636 636 if self.warnings:
637 637 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
638 638 print('maxFreq is setting to %.2f' % (freqrange[-1]))
639 639 maxFreq = freqrange[-1]
640 640
641 641 indminPoint = numpy.where(freqrange >= minFreq)
642 642 indmaxPoint = numpy.where(freqrange <= maxFreq)
643 643
644 644 else:
645 645
646 646 velrange = self.dataOut.getVelRange(1)
647 647
648 648 if minVel == None:
649 649 minVel = velrange[0]
650 650
651 651 if maxVel == None:
652 652 maxVel = velrange[-1]
653 653
654 654 if (minVel < velrange[0]) or (minVel > maxVel):
655 655 if self.warnings:
656 656 print('minVel: %.2f is out of the velocity range' % (minVel))
657 657 print('minVel is setting to %.2f' % (velrange[0]))
658 658 minVel = velrange[0]
659 659
660 660 if (maxVel > velrange[-1]) or (maxVel < minVel):
661 661 if self.warnings:
662 662 print('maxVel: %.2f is out of the velocity range' % (maxVel))
663 663 print('maxVel is setting to %.2f' % (velrange[-1]))
664 664 maxVel = velrange[-1]
665 665
666 666 indminPoint = numpy.where(velrange >= minVel)
667 667 indmaxPoint = numpy.where(velrange <= maxVel)
668 668
669 669
670 670 # seleccion de indices para rango
671 671 minIndex = 0
672 672 maxIndex = 0
673 673 heights = self.dataOut.heightList
674 674
675 675 inda = numpy.where(heights >= minHei)
676 676 indb = numpy.where(heights <= maxHei)
677 677
678 678 try:
679 679 minIndex = inda[0][0]
680 680 except:
681 681 minIndex = 0
682 682
683 683 try:
684 684 maxIndex = indb[0][-1]
685 685 except:
686 686 maxIndex = len(heights)
687 687
688 688 if (minIndex < 0) or (minIndex > maxIndex):
689 689 raise ValueError("some value in (%d,%d) is not valid" % (
690 690 minIndex, maxIndex))
691 691
692 692 if (maxIndex >= self.dataOut.nHeights):
693 693 maxIndex = self.dataOut.nHeights - 1
694 694 #############################################################3
695 695 # seleccion de indices para velocidades
696 696 if self.dataOut.type == 'Spectra':
697 697 try:
698 698 minIndexFFT = indminPoint[0][0]
699 699 except:
700 700 minIndexFFT = 0
701 701
702 702 try:
703 703 maxIndexFFT = indmaxPoint[0][-1]
704 704 except:
705 705 maxIndexFFT = len( self.dataOut.getFreqRange(1))
706 706
707 707 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
708 708 self.isConfig = True
709 709 self.offset = 1
710 710 if offset!=None:
711 711 self.offset = 10**(offset/10)
712 712 #print("config getNoiseB Done")
713 713
714 714 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
715 715 self.dataOut = dataOut
716 716
717 717 if not self.isConfig:
718 718 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
719 719
720 720 self.dataOut.noise_estimation = None
721 721 noise = None
722 722 #print("data type: ",self.dataOut.type, self.dataOut.nIncohInt, self.dataOut.max_nIncohInt)
723 723 if self.dataOut.type == 'Voltage':
724 724 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
725 725 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
726 726 elif self.dataOut.type == 'Spectra':
727 727 #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt)
728 728 noise = numpy.zeros( self.dataOut.nChannels)
729 729 norm = 1
730 730
731 731 for channel in range( self.dataOut.nChannels):
732 732 if not hasattr(self.dataOut.nIncohInt,'__len__'):
733 733 norm = 1
734 734 else:
735 735 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
736 736
737 737 #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape, self.dataOut.max_nIncohInt)
738 738 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
739 739 daux = numpy.multiply(daux, norm)
740 740 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
741 741 # noise[channel] = self.getNoiseByMean(daux)/self.offset
742 742 #print(daux.shape, daux)
743 743 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
744 744 sortdata = numpy.sort(daux, axis=None)
745 745
746 746 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
747 747 #print("noise shape", noise[channel], self.name)
748 748
749 749 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
750 750 else:
751 751 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
752 752
753 753 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
754 754 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
755 755 #print("2: ",self.dataOut.noise_estimation)
756 756 #print(self.dataOut.flagNoData)
757 757 #print("getNoise Done", 10*numpy.log10(noise))
758 758 return self.dataOut
759 759
760 760 def getNoiseByMean(self,data):
761 761 #data debe estar ordenado
762 762 data = numpy.mean(data,axis=1)
763 763 sortdata = numpy.sort(data, axis=None)
764 764 #sortID=data.argsort()
765 765 #print(data.shape)
766 766
767 767 pnoise = None
768 768 j = 0
769 769
770 770 mean = numpy.mean(sortdata)
771 771 min = numpy.min(sortdata)
772 772 delta = mean - min
773 773 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
774 774 #print(len(indexes))
775 775 if len(indexes)==0:
776 776 pnoise = numpy.mean(sortdata)
777 777 else:
778 778 j = indexes[0]
779 779 pnoise = numpy.mean(sortdata[0:j])
780 780
781 781 # from matplotlib import pyplot as plt
782 782 # plt.plot(sortdata)
783 783 # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r')
784 784 # plt.show()
785 785 #print("noise: ", 10*numpy.log10(pnoise))
786 786 return pnoise
787 787
788 788 def getNoiseByHS(self,data, navg):
789 789 #data debe estar ordenado
790 790 #data = numpy.mean(data,axis=1)
791 791 sortdata = numpy.sort(data, axis=None)
792 792
793 793 lenOfData = len(sortdata)
794 794 nums_min = lenOfData*0.2
795 795
796 796 if nums_min <= 5:
797 797
798 798 nums_min = 5
799 799
800 800 sump = 0.
801 801 sumq = 0.
802 802
803 803 j = 0
804 804 cont = 1
805 805
806 806 while((cont == 1)and(j < lenOfData)):
807 807
808 808 sump += sortdata[j]
809 809 sumq += sortdata[j]**2
810 810 #sumq -= sump**2
811 811 if j > nums_min:
812 812 rtest = float(j)/(j-1) + 1.0/navg
813 813 #if ((sumq*j) > (sump**2)):
814 814 if ((sumq*j) > (rtest*sump**2)):
815 815 j = j - 1
816 816 sump = sump - sortdata[j]
817 817 sumq = sumq - sortdata[j]**2
818 818 cont = 0
819 819
820 820 j += 1
821 821
822 822 lnoise = sump / j
823 823
824 824 return lnoise
825 825
826 826
827 827
828 828 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
829 829 z = (x - a1) / a2
830 830 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
831 831 return y
832 832
833 833
834 834 # class CleanRayleigh(Operation):
835 835 #
836 836 # def __init__(self):
837 837 #
838 838 # Operation.__init__(self)
839 839 # self.i=0
840 840 # self.isConfig = False
841 841 # self.__dataReady = False
842 842 # self.__profIndex = 0
843 843 # self.byTime = False
844 844 # self.byProfiles = False
845 845 #
846 846 # self.bloques = None
847 847 # self.bloque0 = None
848 848 #
849 849 # self.index = 0
850 850 #
851 851 # self.buffer = 0
852 852 # self.buffer2 = 0
853 853 # self.buffer3 = 0
854 854 #
855 855 #
856 856 # def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
857 857 #
858 858 # self.nChannels = dataOut.nChannels
859 859 # self.nProf = dataOut.nProfiles
860 860 # self.nPairs = dataOut.data_cspc.shape[0]
861 861 # self.pairsArray = numpy.array(dataOut.pairsList)
862 862 # self.spectra = dataOut.data_spc
863 863 # self.cspectra = dataOut.data_cspc
864 864 # self.heights = dataOut.heightList #alturas totales
865 865 # self.nHeights = len(self.heights)
866 866 # self.min_hei = min_hei
867 867 # self.max_hei = max_hei
868 868 # if (self.min_hei == None):
869 869 # self.min_hei = 0
870 870 # if (self.max_hei == None):
871 871 # self.max_hei = dataOut.heightList[-1]
872 872 # self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
873 873 # self.heightsClean = self.heights[self.hval] #alturas filtradas
874 874 # self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
875 875 # self.nHeightsClean = len(self.heightsClean)
876 876 # self.channels = dataOut.channelList
877 877 # self.nChan = len(self.channels)
878 878 # self.nIncohInt = dataOut.nIncohInt
879 879 # self.__initime = dataOut.utctime
880 880 # self.maxAltInd = self.hval[-1]+1
881 881 # self.minAltInd = self.hval[0]
882 882 #
883 883 # self.crosspairs = dataOut.pairsList
884 884 # self.nPairs = len(self.crosspairs)
885 885 # self.normFactor = dataOut.normFactor
886 886 # self.nFFTPoints = dataOut.nFFTPoints
887 887 # self.ippSeconds = dataOut.ippSeconds
888 888 # self.currentTime = self.__initime
889 889 # self.pairsArray = numpy.array(dataOut.pairsList)
890 890 # self.factor_stdv = factor_stdv
891 891 #
892 892 # if n != None :
893 893 # self.byProfiles = True
894 894 # self.nIntProfiles = n
895 895 # else:
896 896 # self.__integrationtime = timeInterval
897 897 #
898 898 # self.__dataReady = False
899 899 # self.isConfig = True
900 900 #
901 901 #
902 902 #
903 903 # def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
904 904 # #print("runing cleanRayleigh")
905 905 # if not self.isConfig :
906 906 #
907 907 # self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
908 908 #
909 909 # tini=dataOut.utctime
910 910 #
911 911 # if self.byProfiles:
912 912 # if self.__profIndex == self.nIntProfiles:
913 913 # self.__dataReady = True
914 914 # else:
915 915 # if (tini - self.__initime) >= self.__integrationtime:
916 916 #
917 917 # self.__dataReady = True
918 918 # self.__initime = tini
919 919 #
920 920 # #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
921 921 #
922 922 # if self.__dataReady:
923 923 #
924 924 # self.__profIndex = 0
925 925 # jspc = self.buffer
926 926 # jcspc = self.buffer2
927 927 # #jnoise = self.buffer3
928 928 # self.buffer = dataOut.data_spc
929 929 # self.buffer2 = dataOut.data_cspc
930 930 # #self.buffer3 = dataOut.noise
931 931 # self.currentTime = dataOut.utctime
932 932 # if numpy.any(jspc) :
933 933 # #print( jspc.shape, jcspc.shape)
934 934 # jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
935 935 # try:
936 936 # jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
937 937 # except:
938 938 # print("no cspc")
939 939 # self.__dataReady = False
940 940 # #print( jspc.shape, jcspc.shape)
941 941 # dataOut.flagNoData = False
942 942 # else:
943 943 # dataOut.flagNoData = True
944 944 # self.__dataReady = False
945 945 # return dataOut
946 946 # else:
947 947 # #print( len(self.buffer))
948 948 # if numpy.any(self.buffer):
949 949 # self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
950 950 # try:
951 951 # self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
952 952 # self.buffer3 += dataOut.data_dc
953 953 # except:
954 954 # pass
955 955 # else:
956 956 # self.buffer = dataOut.data_spc
957 957 # self.buffer2 = dataOut.data_cspc
958 958 # self.buffer3 = dataOut.data_dc
959 959 # #print self.index, self.fint
960 960 # #print self.buffer2.shape
961 961 # dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
962 962 # self.__profIndex += 1
963 963 # return dataOut ## NOTE: REV
964 964 #
965 965 #
966 966 # #index = tini.tm_hour*12+tini.tm_min/5
967 967 # '''
968 968 # #REVISAR
969 969 # '''
970 970 # # jspc = jspc/self.nFFTPoints/self.normFactor
971 971 # # jcspc = jcspc/self.nFFTPoints/self.normFactor
972 972 #
973 973 #
974 974 #
975 975 # tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
976 976 # dataOut.data_spc = tmp_spectra
977 977 # dataOut.data_cspc = tmp_cspectra
978 978 #
979 979 # #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
980 980 #
981 981 # dataOut.data_dc = self.buffer3
982 982 # dataOut.nIncohInt *= self.nIntProfiles
983 983 # dataOut.max_nIncohInt = self.nIntProfiles
984 984 # dataOut.utctime = self.currentTime #tiempo promediado
985 985 # #print("Time: ",time.localtime(dataOut.utctime))
986 986 # # dataOut.data_spc = sat_spectra
987 987 # # dataOut.data_cspc = sat_cspectra
988 988 # self.buffer = 0
989 989 # self.buffer2 = 0
990 990 # self.buffer3 = 0
991 991 #
992 992 # return dataOut
993 993 #
994 994 # def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
995 995 # print("OP cleanRayleigh")
996 996 # #import matplotlib.pyplot as plt
997 997 # #for k in range(149):
998 998 # #channelsProcssd = []
999 999 # #channelA_ok = False
1000 1000 # #rfunc = cspectra.copy() #self.bloques
1001 1001 # rfunc = spectra.copy()
1002 1002 # #rfunc = cspectra
1003 1003 # #val_spc = spectra*0.0 #self.bloque0*0.0
1004 1004 # #val_cspc = cspectra*0.0 #self.bloques*0.0
1005 1005 # #in_sat_spectra = spectra.copy() #self.bloque0
1006 1006 # #in_sat_cspectra = cspectra.copy() #self.bloques
1007 1007 #
1008 1008 #
1009 1009 # ###ONLY FOR TEST:
1010 1010 # raxs = math.ceil(math.sqrt(self.nPairs))
1011 1011 # if raxs == 0:
1012 1012 # raxs = 1
1013 1013 # caxs = math.ceil(self.nPairs/raxs)
1014 1014 # if self.nPairs <4:
1015 1015 # raxs = 2
1016 1016 # caxs = 2
1017 1017 # #print(raxs, caxs)
1018 1018 # fft_rev = 14 #nFFT to plot
1019 1019 # hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
1020 1020 # hei_rev = hei_rev[0]
1021 1021 # #print(hei_rev)
1022 1022 #
1023 1023 # #print numpy.absolute(rfunc[:,0,0,14])
1024 1024 #
1025 1025 # gauss_fit, covariance = None, None
1026 1026 # for ih in range(self.minAltInd,self.maxAltInd):
1027 1027 # for ifreq in range(self.nFFTPoints):
1028 1028 # '''
1029 1029 # ###ONLY FOR TEST:
1030 1030 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1031 1031 # fig, axs = plt.subplots(raxs, caxs)
1032 1032 # fig2, axs2 = plt.subplots(raxs, caxs)
1033 1033 # col_ax = 0
1034 1034 # row_ax = 0
1035 1035 # '''
1036 1036 # #print(self.nPairs)
1037 1037 # for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
1038 1038 # # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
1039 1039 # # continue
1040 1040 # # if not self.crosspairs[ii][0] in channelsProcssd:
1041 1041 # # channelA_ok = True
1042 1042 # #print("pair: ",self.crosspairs[ii])
1043 1043 # '''
1044 1044 # ###ONLY FOR TEST:
1045 1045 # if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
1046 1046 # col_ax = 0
1047 1047 # row_ax += 1
1048 1048 # '''
1049 1049 # func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
1050 1050 # #print(func2clean.shape)
1051 1051 # val = (numpy.isfinite(func2clean)==True).nonzero()
1052 1052 #
1053 1053 # if len(val)>0: #limitador
1054 1054 # min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
1055 1055 # if min_val <= -40 :
1056 1056 # min_val = -40
1057 1057 # max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
1058 1058 # if max_val >= 200 :
1059 1059 # max_val = 200
1060 1060 # #print min_val, max_val
1061 1061 # step = 1
1062 1062 # #print("Getting bins and the histogram")
1063 1063 # x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
1064 1064 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1065 1065 # #print(len(y_dist),len(binstep[:-1]))
1066 1066 # #print(row_ax,col_ax, " ..")
1067 1067 # #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
1068 1068 # mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
1069 1069 # sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
1070 1070 # parg = [numpy.amax(y_dist),mean,sigma]
1071 1071 #
1072 1072 # newY = None
1073 1073 #
1074 1074 # try :
1075 1075 # gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
1076 1076 # mode = gauss_fit[1]
1077 1077 # stdv = gauss_fit[2]
1078 1078 # #print(" FIT OK",gauss_fit)
1079 1079 # '''
1080 1080 # ###ONLY FOR TEST:
1081 1081 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1082 1082 # newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
1083 1083 # axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1084 1084 # axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1085 1085 # axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1086 1086 # '''
1087 1087 # except:
1088 1088 # mode = mean
1089 1089 # stdv = sigma
1090 1090 # #print("FIT FAIL")
1091 1091 # #continue
1092 1092 #
1093 1093 #
1094 1094 # #print(mode,stdv)
1095 1095 # #Removing echoes greater than mode + std_factor*stdv
1096 1096 # noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
1097 1097 # #noval tiene los indices que se van a remover
1098 1098 # #print("Chan ",ii," novals: ",len(noval[0]))
1099 1099 # if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
1100 1100 # novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
1101 1101 # #print(novall)
1102 1102 # #print(" ",self.pairsArray[ii])
1103 1103 # #cross_pairs = self.pairsArray[ii]
1104 1104 # #Getting coherent echoes which are removed.
1105 1105 # # if len(novall[0]) > 0:
1106 1106 # #
1107 1107 # # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
1108 1108 # # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
1109 1109 # # val_cspc[novall[0],ii,ifreq,ih] = 1
1110 1110 # #print("OUT NOVALL 1")
1111 1111 # try:
1112 1112 # pair = (self.channels[ii],self.channels[ii + 1])
1113 1113 # except:
1114 1114 # pair = (99,99)
1115 1115 # #print("par ", pair)
1116 1116 # if ( pair in self.crosspairs):
1117 1117 # q = self.crosspairs.index(pair)
1118 1118 # #print("estΓ‘ aqui: ", q, (ii,ii + 1))
1119 1119 # new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
1120 1120 # cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
1121 1121 #
1122 1122 # #if channelA_ok:
1123 1123 # #chA = self.channels.index(cross_pairs[0])
1124 1124 # new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
1125 1125 # spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
1126 1126 # #channelA_ok = False
1127 1127 #
1128 1128 # # chB = self.channels.index(cross_pairs[1])
1129 1129 # # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
1130 1130 # # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
1131 1131 # #
1132 1132 # # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
1133 1133 # # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
1134 1134 # '''
1135 1135 # ###ONLY FOR TEST:
1136 1136 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1137 1137 # func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
1138 1138 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1139 1139 # axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1140 1140 # axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1141 1141 # axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1142 1142 # '''
1143 1143 # '''
1144 1144 # ###ONLY FOR TEST:
1145 1145 # col_ax += 1 #contador de ploteo columnas
1146 1146 # ##print(col_ax)
1147 1147 # ###ONLY FOR TEST:
1148 1148 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1149 1149 # title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
1150 1150 # title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
1151 1151 # fig.suptitle(title)
1152 1152 # fig2.suptitle(title2)
1153 1153 # plt.show()
1154 1154 # '''
1155 1155 # ##################################################################################################
1156 1156 #
1157 1157 # #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
1158 1158 # out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
1159 1159 # out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
1160 1160 # for ih in range(self.nHeights):
1161 1161 # for ifreq in range(self.nFFTPoints):
1162 1162 # for ich in range(self.nChan):
1163 1163 # tmp = spectra[:,ich,ifreq,ih]
1164 1164 # valid = (numpy.isfinite(tmp[:])==True).nonzero()
1165 1165 #
1166 1166 # if len(valid[0]) >0 :
1167 1167 # out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1168 1168 #
1169 1169 # for icr in range(self.nPairs):
1170 1170 # tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
1171 1171 # valid = (numpy.isfinite(tmp)==True).nonzero()
1172 1172 # if len(valid[0]) > 0:
1173 1173 # out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1174 1174 #
1175 1175 # return out_spectra, out_cspectra
1176 1176 #
1177 1177 # def REM_ISOLATED_POINTS(self,array,rth):
1178 1178 # # import matplotlib.pyplot as plt
1179 1179 # if rth == None :
1180 1180 # rth = 4
1181 1181 # #print("REM ISO")
1182 1182 # num_prof = len(array[0,:,0])
1183 1183 # num_hei = len(array[0,0,:])
1184 1184 # n2d = len(array[:,0,0])
1185 1185 #
1186 1186 # for ii in range(n2d) :
1187 1187 # #print ii,n2d
1188 1188 # tmp = array[ii,:,:]
1189 1189 # #print tmp.shape, array[ii,101,:],array[ii,102,:]
1190 1190 #
1191 1191 # # fig = plt.figure(figsize=(6,5))
1192 1192 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1193 1193 # # ax = fig.add_axes([left, bottom, width, height])
1194 1194 # # x = range(num_prof)
1195 1195 # # y = range(num_hei)
1196 1196 # # cp = ax.contour(y,x,tmp)
1197 1197 # # ax.clabel(cp, inline=True,fontsize=10)
1198 1198 # # plt.show()
1199 1199 #
1200 1200 # #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
1201 1201 # tmp = numpy.reshape(tmp,num_prof*num_hei)
1202 1202 # indxs1 = (numpy.isfinite(tmp)==True).nonzero()
1203 1203 # indxs2 = (tmp > 0).nonzero()
1204 1204 #
1205 1205 # indxs1 = (indxs1[0])
1206 1206 # indxs2 = indxs2[0]
1207 1207 # #indxs1 = numpy.array(indxs1[0])
1208 1208 # #indxs2 = numpy.array(indxs2[0])
1209 1209 # indxs = None
1210 1210 # #print indxs1 , indxs2
1211 1211 # for iv in range(len(indxs2)):
1212 1212 # indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
1213 1213 # #print len(indxs2), indv
1214 1214 # if len(indv[0]) > 0 :
1215 1215 # indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
1216 1216 # # print indxs
1217 1217 # indxs = indxs[1:]
1218 1218 # #print(indxs, len(indxs))
1219 1219 # if len(indxs) < 4 :
1220 1220 # array[ii,:,:] = 0.
1221 1221 # return
1222 1222 #
1223 1223 # xpos = numpy.mod(indxs ,num_hei)
1224 1224 # ypos = (indxs / num_hei)
1225 1225 # sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1226 1226 # #print sx
1227 1227 # xpos = xpos[sx]
1228 1228 # ypos = ypos[sx]
1229 1229 #
1230 1230 # # *********************************** Cleaning isolated points **********************************
1231 1231 # ic = 0
1232 1232 # while True :
1233 1233 # r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1234 1234 # #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1235 1235 # #plt.plot(r)
1236 1236 # #plt.show()
1237 1237 # no_coh1 = (numpy.isfinite(r)==True).nonzero()
1238 1238 # no_coh2 = (r <= rth).nonzero()
1239 1239 # #print r, no_coh1, no_coh2
1240 1240 # no_coh1 = numpy.array(no_coh1[0])
1241 1241 # no_coh2 = numpy.array(no_coh2[0])
1242 1242 # no_coh = None
1243 1243 # #print valid1 , valid2
1244 1244 # for iv in range(len(no_coh2)):
1245 1245 # indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1246 1246 # if len(indv[0]) > 0 :
1247 1247 # no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1248 1248 # no_coh = no_coh[1:]
1249 1249 # #print len(no_coh), no_coh
1250 1250 # if len(no_coh) < 4 :
1251 1251 # #print xpos[ic], ypos[ic], ic
1252 1252 # # plt.plot(r)
1253 1253 # # plt.show()
1254 1254 # xpos[ic] = numpy.nan
1255 1255 # ypos[ic] = numpy.nan
1256 1256 #
1257 1257 # ic = ic + 1
1258 1258 # if (ic == len(indxs)) :
1259 1259 # break
1260 1260 # #print( xpos, ypos)
1261 1261 #
1262 1262 # indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1263 1263 # #print indxs[0]
1264 1264 # if len(indxs[0]) < 4 :
1265 1265 # array[ii,:,:] = 0.
1266 1266 # return
1267 1267 #
1268 1268 # xpos = xpos[indxs[0]]
1269 1269 # ypos = ypos[indxs[0]]
1270 1270 # for i in range(0,len(ypos)):
1271 1271 # ypos[i]=int(ypos[i])
1272 1272 # junk = tmp
1273 1273 # tmp = junk*0.0
1274 1274 #
1275 1275 # tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1276 1276 # array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1277 1277 #
1278 1278 # #print array.shape
1279 1279 # #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1280 1280 # #print tmp.shape
1281 1281 #
1282 1282 # # fig = plt.figure(figsize=(6,5))
1283 1283 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1284 1284 # # ax = fig.add_axes([left, bottom, width, height])
1285 1285 # # x = range(num_prof)
1286 1286 # # y = range(num_hei)
1287 1287 # # cp = ax.contour(y,x,array[ii,:,:])
1288 1288 # # ax.clabel(cp, inline=True,fontsize=10)
1289 1289 # # plt.show()
1290 1290 # return array
1291 1291 #
1292 1292
1293 1293 class IntegrationFaradaySpectra(Operation):
1294 1294
1295 1295 __profIndex = 0
1296 1296 __withOverapping = False
1297 1297
1298 1298 __byTime = False
1299 1299 __initime = None
1300 1300 __lastdatatime = None
1301 1301 __integrationtime = None
1302 1302
1303 1303 __buffer_spc = None
1304 1304 __buffer_cspc = None
1305 1305 __buffer_dc = None
1306 1306
1307 1307 __dataReady = False
1308 1308
1309 1309 __timeInterval = None
1310 1310 n_ints = None #matriz de numero de integracions (CH,HEI)
1311 1311 n = None
1312 1312 minHei_ind = None
1313 1313 maxHei_ind = None
1314 1314 navg = 1.0
1315 1315 factor = 0.0
1316 1316 dataoutliers = None # (CHANNELS, HEIGHTS)
1317 1317
1318 1318 _flagProfilesByRange = False
1319 1319 _nProfilesByRange = 0
1320 1320
1321 1321 def __init__(self):
1322 1322
1323 1323 Operation.__init__(self)
1324 1324
1325 1325 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1326 1326 """
1327 1327 Set the parameters of the integration class.
1328 1328
1329 1329 Inputs:
1330 1330
1331 1331 n : Number of coherent integrations
1332 1332 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1333 1333 overlapping :
1334 1334
1335 1335 """
1336 1336
1337 1337 self.__initime = None
1338 1338 self.__lastdatatime = 0
1339 1339
1340 1340 self.__buffer_spc = []
1341 1341 self.__buffer_cspc = []
1342 1342 self.__buffer_dc = 0
1343 1343
1344 1344 self.__profIndex = 0
1345 1345 self.__dataReady = False
1346 1346 self.__byTime = False
1347 1347
1348 1348 self.factor = factor
1349 1349 self.navg = avg
1350 1350 #self.ByLags = dataOut.ByLags ###REDEFINIR
1351 1351 self.ByLags = False
1352 1352 self.maxProfilesInt = 0
1353 1353 self.__nChannels = dataOut.nChannels
1354 1354 if DPL != None:
1355 1355 self.DPL=DPL
1356 1356 else:
1357 1357 #self.DPL=dataOut.DPL ###REDEFINIR
1358 1358 self.DPL=0
1359 1359
1360 1360 if n is None and timeInterval is None:
1361 1361 raise ValueError("n or timeInterval should be specified ...")
1362 1362
1363 1363 if n is not None:
1364 1364 self.n = int(n)
1365 1365 else:
1366 1366 self.__integrationtime = int(timeInterval)
1367 1367 self.n = None
1368 1368 self.__byTime = True
1369 1369
1370 1370
1371 1371 if minHei == None:
1372 1372 minHei = self.dataOut.heightList[0]
1373 1373
1374 1374 if maxHei == None:
1375 1375 maxHei = self.dataOut.heightList[-1]
1376 1376
1377 1377 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
1378 1378 print('minHei: %.2f is out of the heights range' % (minHei))
1379 1379 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
1380 1380 minHei = self.dataOut.heightList[0]
1381 1381
1382 1382 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
1383 1383 print('maxHei: %.2f is out of the heights range' % (maxHei))
1384 1384 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
1385 1385 maxHei = self.dataOut.heightList[-1]
1386 1386
1387 1387 ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
1388 1388 ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
1389 1389 self.minHei_ind = ind_list1[0][0]
1390 1390 self.maxHei_ind = ind_list2[0][-1]
1391 1391
1392 1392 def putData(self, data_spc, data_cspc, data_dc):
1393 1393 """
1394 1394 Add a profile to the __buffer_spc and increase in one the __profileIndex
1395 1395
1396 1396 """
1397 1397
1398 1398 self.__buffer_spc.append(data_spc)
1399 1399
1400 1400 if self.__nChannels < 2:
1401 1401 self.__buffer_cspc = None
1402 1402 else:
1403 1403 self.__buffer_cspc.append(data_cspc)
1404 1404
1405 1405 if data_dc is None:
1406 1406 self.__buffer_dc = None
1407 1407 else:
1408 1408 self.__buffer_dc += data_dc
1409 1409
1410 1410 self.__profIndex += 1
1411 1411
1412 1412 return
1413 1413
1414 1414 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1415 1415 #data debe estar ordenado
1416 1416 #sortdata = numpy.sort(data, axis=None)
1417 1417 #sortID=data.argsort()
1418 1418 lenOfData = len(sortdata)
1419 1419 nums_min = lenOfData*factor
1420 1420 if nums_min <= 5:
1421 1421 nums_min = 5
1422 1422 sump = 0.
1423 1423 sumq = 0.
1424 1424 j = 0
1425 1425 cont = 1
1426 1426 while((cont == 1)and(j < lenOfData)):
1427 1427 sump += sortdata[j]
1428 1428 sumq += sortdata[j]**2
1429 1429 if j > nums_min:
1430 1430 rtest = float(j)/(j-1) + 1.0/navg
1431 1431 if ((sumq*j) > (rtest*sump**2)):
1432 1432 j = j - 1
1433 1433 sump = sump - sortdata[j]
1434 1434 sumq = sumq - sortdata[j]**2
1435 1435 cont = 0
1436 1436 j += 1
1437 1437 #lnoise = sump / j
1438 1438 #print("H S done")
1439 1439 #return j,sortID
1440 1440 return j
1441 1441
1442 1442
1443 1443 def pushData(self):
1444 1444 """
1445 1445 Return the sum of the last profiles and the profiles used in the sum.
1446 1446
1447 1447 Affected:
1448 1448
1449 1449 self.__profileIndex
1450 1450
1451 1451 """
1452 1452 bufferH=None
1453 1453 buffer=None
1454 1454 buffer1=None
1455 1455 buffer_cspc=None
1456 1456 #print("aes: ", self.__buffer_cspc)
1457 1457 self.__buffer_spc=numpy.array(self.__buffer_spc)
1458 1458 if self.__nChannels > 1 :
1459 1459 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1460 1460
1461 1461 #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
1462 1462
1463 1463 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1464 1464 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1465 1465
1466 1466 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1467 1467
1468 1468 for k in range(self.minHei_ind,self.maxHei_ind):
1469 1469 if self.__nChannels > 1:
1470 1470 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1471 1471
1472 1472 outliers_IDs_cspc=[]
1473 1473 cspc_outliers_exist=False
1474 1474 for i in range(self.nChannels):#dataOut.nChannels):
1475 1475
1476 1476 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1477 1477 indexes=[]
1478 1478 #sortIDs=[]
1479 1479 outliers_IDs=[]
1480 1480
1481 1481 for j in range(self.nProfiles): #frecuencias en el tiempo
1482 1482 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1483 1483 # continue
1484 1484 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1485 1485 # continue
1486 1486 buffer=buffer1[:,j]
1487 1487 sortdata = numpy.sort(buffer, axis=None)
1488 1488
1489 1489 sortID=buffer.argsort()
1490 1490 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1491 1491
1492 1492 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1493 1493
1494 1494 # fig,ax = plt.subplots()
1495 1495 # ax.set_title(str(k)+" "+str(j))
1496 1496 # x=range(len(sortdata))
1497 1497 # ax.scatter(x,sortdata)
1498 1498 # ax.axvline(index)
1499 1499 # plt.show()
1500 1500
1501 1501 indexes.append(index)
1502 1502 #sortIDs.append(sortID)
1503 1503 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1504 1504
1505 1505 #print("Outliers: ",outliers_IDs)
1506 1506 outliers_IDs=numpy.array(outliers_IDs)
1507 1507 outliers_IDs=outliers_IDs.ravel()
1508 1508 outliers_IDs=numpy.unique(outliers_IDs)
1509 1509 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1510 1510 indexes=numpy.array(indexes)
1511 1511 indexmin=numpy.min(indexes)
1512 1512
1513 1513
1514 1514 #print(indexmin,buffer1.shape[0], k)
1515 1515
1516 1516 # fig,ax = plt.subplots()
1517 1517 # ax.plot(sortdata)
1518 1518 # ax2 = ax.twinx()
1519 1519 # x=range(len(indexes))
1520 1520 # #plt.scatter(x,indexes)
1521 1521 # ax2.scatter(x,indexes)
1522 1522 # plt.show()
1523 1523
1524 1524 if indexmin != buffer1.shape[0]:
1525 1525 if self.__nChannels > 1:
1526 1526 cspc_outliers_exist= True
1527 1527
1528 1528 lt=outliers_IDs
1529 1529 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1530 1530
1531 1531 for p in list(outliers_IDs):
1532 1532 #buffer1[p,:]=avg
1533 1533 buffer1[p,:] = numpy.NaN
1534 1534
1535 1535 self.dataOutliers[i,k] = len(outliers_IDs)
1536 1536
1537 1537
1538 1538 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1539 1539
1540 1540
1541 1541 if self.__nChannels > 1:
1542 1542 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1543 1543
1544 1544
1545 1545 if self.__nChannels > 1:
1546 1546 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1547 1547 if cspc_outliers_exist:
1548 1548
1549 1549 lt=outliers_IDs_cspc
1550 1550
1551 1551 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1552 1552 for p in list(outliers_IDs_cspc):
1553 1553 #buffer_cspc[p,:]=avg
1554 1554 buffer_cspc[p,:] = numpy.NaN
1555 1555
1556 1556 if self.__nChannels > 1:
1557 1557 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1558 1558
1559 1559
1560 1560
1561 1561
1562 1562 nOutliers = len(outliers_IDs)
1563 1563 #print("Outliers n: ",self.dataOutliers,nOutliers)
1564 1564 buffer=None
1565 1565 bufferH=None
1566 1566 buffer1=None
1567 1567 buffer_cspc=None
1568 1568
1569 1569
1570 1570 buffer=None
1571 1571
1572 1572 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1573 1573 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1574 1574 if self.__nChannels > 1:
1575 1575 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1576 1576 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1577 1577 else:
1578 1578 data_cspc = None
1579 1579 data_dc = self.__buffer_dc
1580 1580 #(CH, HEIGH)
1581 1581 self.maxProfilesInt = self.__profIndex - 1
1582 1582 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1583 1583
1584 1584 self.__buffer_spc = []
1585 1585 self.__buffer_cspc = []
1586 1586 self.__buffer_dc = 0
1587 1587 self.__profIndex = 0
1588 1588 #print("cleaned ",data_cspc)
1589 1589 return data_spc, data_cspc, data_dc, n
1590 1590
1591 1591 def byProfiles(self, *args):
1592 1592
1593 1593 self.__dataReady = False
1594 1594 avgdata_spc = None
1595 1595 avgdata_cspc = None
1596 1596 avgdata_dc = None
1597 1597
1598 1598 self.putData(*args)
1599 1599
1600 1600 if self.__profIndex >= self.n:
1601 1601
1602 1602 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1603 1603 self.n_ints = n
1604 1604 self.__dataReady = True
1605 1605
1606 1606 return avgdata_spc, avgdata_cspc, avgdata_dc
1607 1607
1608 1608 def byTime(self, datatime, *args):
1609 1609
1610 1610 self.__dataReady = False
1611 1611 avgdata_spc = None
1612 1612 avgdata_cspc = None
1613 1613 avgdata_dc = None
1614 1614
1615 1615 self.putData(*args)
1616 1616
1617 1617 if (datatime - self.__initime) >= self.__integrationtime:
1618 1618 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1619 1619 self.n_ints = n
1620 1620 self.__dataReady = True
1621 1621
1622 1622 return avgdata_spc, avgdata_cspc, avgdata_dc
1623 1623
1624 1624 def integrate(self, datatime, *args):
1625 1625
1626 1626 if self.__profIndex == 0:
1627 1627 self.__initime = datatime
1628 1628
1629 1629 if self.__byTime:
1630 1630 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1631 1631 datatime, *args)
1632 1632 else:
1633 1633 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1634 1634
1635 1635 if not self.__dataReady:
1636 1636 return None, None, None, None
1637 1637
1638 1638 #print("integrate", avgdata_cspc)
1639 1639 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1640 1640
1641 1641 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1642 1642 self.dataOut = dataOut
1643 1643 if n == 1:
1644 1644 return self.dataOut
1645 1645 self.dataOut.processingHeaderObj.timeIncohInt = timeInterval
1646 1646
1647 1647 if dataOut.flagProfilesByRange:
1648 1648 self._flagProfilesByRange = True
1649 1649
1650 1650 if self.dataOut.nChannels == 1:
1651 1651 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1652 1652 #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1653 1653 if not self.isConfig:
1654 1654 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1655 1655 self.isConfig = True
1656 1656
1657 1657 if not self.ByLags:
1658 1658 self.nProfiles=self.dataOut.nProfiles
1659 1659 self.nChannels=self.dataOut.nChannels
1660 1660 self.nHeights=self.dataOut.nHeights
1661 1661 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1662 1662 self.dataOut.data_spc,
1663 1663 self.dataOut.data_cspc,
1664 1664 self.dataOut.data_dc)
1665 1665 else:
1666 1666 self.nProfiles=self.dataOut.nProfiles
1667 1667 self.nChannels=self.dataOut.nChannels
1668 1668 self.nHeights=self.dataOut.nHeights
1669 1669 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
1670 1670 self.dataOut.dataLag_spc,
1671 1671 self.dataOut.dataLag_cspc,
1672 1672 self.dataOut.dataLag_dc)
1673 1673 self.dataOut.flagNoData = True
1674 1674
1675 1675 if self._flagProfilesByRange:
1676 1676 dataOut.flagProfilesByRange = True
1677 1677 self._nProfilesByRange += dataOut.nProfilesByRange
1678 1678
1679 1679 if self.__dataReady:
1680 1680
1681 1681 if not self.ByLags:
1682 1682 if self.nChannels == 1:
1683 1683 #print("f int", avgdata_spc.shape)
1684 1684 self.dataOut.data_spc = avgdata_spc
1685 1685 self.dataOut.data_cspc = None
1686 1686 else:
1687 1687 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1688 1688 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1689 1689 self.dataOut.data_dc = avgdata_dc
1690 1690 self.dataOut.data_outlier = self.dataOutliers
1691 1691
1692 1692
1693 1693 else:
1694 1694 self.dataOut.dataLag_spc = avgdata_spc
1695 1695 self.dataOut.dataLag_cspc = avgdata_cspc
1696 1696 self.dataOut.dataLag_dc = avgdata_dc
1697 1697
1698 1698 self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
1699 1699 self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
1700 1700 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1701 1701
1702 1702 self.dataOut.nIncohInt *= self.n_ints
1703 1703 #print("maxProfilesInt: ",self.maxProfilesInt)
1704 1704
1705 1705 self.dataOut.utctime = avgdatatime
1706 1706 self.dataOut.flagNoData = False
1707 1707
1708 1708 dataOut.nProfilesByRange = self._nProfilesByRange
1709 1709 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
1710 1710 self._flagProfilesByRange = False
1711 1711
1712 1712 # #update Processing Header:
1713 1713 # self.dataOut.processingHeaderObj.nIncohInt =
1714 1714 # self.dataOut.processingHeaderObj.nFFTPoints = self.dataOut.nFFTPoints
1715 1715
1716 1716 #print("Faraday Integration DONE...", self.dataOut.data_cspc)
1717 1717 #print(self.dataOut.flagNoData)
1718 1718 return self.dataOut
1719 1719
1720 1720
1721 1721
1722 1722 class removeInterference(Operation):
1723 1723
1724 1724 def removeInterference3(self, min_hei = None, max_hei = None):
1725 1725
1726 1726 jspectra = self.dataOut.data_spc
1727 1727 #jcspectra = self.dataOut.data_cspc
1728 1728 jnoise = self.dataOut.getNoise()
1729 1729 num_incoh = self.dataOut.max_nIncohInt
1730 1730 #print(jspectra.shape)
1731 1731 num_channel, num_prof, num_hei = jspectra.shape
1732 1732 minHei = min_hei
1733 1733 maxHei = max_hei
1734 1734 ########################################################################
1735 1735 if minHei == None or (minHei < self.dataOut.heightList[0]):
1736 1736 minHei = self.dataOut.heightList[0]
1737 1737
1738 1738 if maxHei == None or (maxHei > self.dataOut.heightList[-1]):
1739 1739 maxHei = self.dataOut.heightList[-1]
1740 1740 minIndex = 0
1741 1741 maxIndex = 0
1742 1742 heights = self.dataOut.heightList
1743 1743
1744 1744 inda = numpy.where(heights >= minHei)
1745 1745 indb = numpy.where(heights <= maxHei)
1746 1746
1747 1747 try:
1748 1748 minIndex = inda[0][0]
1749 1749 except:
1750 1750 minIndex = 0
1751 1751 try:
1752 1752 maxIndex = indb[0][-1]
1753 1753 except:
1754 1754 maxIndex = len(heights)
1755 1755
1756 1756 if (minIndex < 0) or (minIndex > maxIndex):
1757 1757 raise ValueError("some value in (%d,%d) is not valid" % (
1758 1758 minIndex, maxIndex))
1759 1759 if (maxIndex >= self.dataOut.nHeights):
1760 1760 maxIndex = self.dataOut.nHeights - 1
1761 1761
1762 1762 ########################################################################
1763 1763
1764 1764
1765 1765 #dataOut.max_nIncohInt * dataOut.nCohInt
1766 1766 if hasattr(self.dataOut.nIncohInt, 'shape'):
1767 1767 norm = self.dataOut.nIncohInt.T /self.dataOut.max_nIncohInt
1768 1768 norm = norm.T
1769 1769 else:
1770 1770 norm = self.dataOut.nIncohInt /self.dataOut.max_nIncohInt
1771 1771 norm = norm
1772 1772
1773 1773 # Subrutina de Remocion de la Interferencia
1774 1774 for ich in range(num_channel):
1775 1775 # Se ordena los espectros segun su potencia (menor a mayor)
1776 1776 #power = jspectra[ich, mask_prof, :]
1777 1777 if hasattr(self.dataOut.nIncohInt, 'shape'):
1778 1778 interf = jspectra[ich, :, minIndex:maxIndex]/norm[ich,minIndex:maxIndex]
1779 1779 else:
1780 1780 interf = jspectra[ich, :, minIndex:maxIndex]/norm
1781 1781 # print(interf.shape)
1782 1782 inttef = interf.mean(axis=1)
1783 1783
1784 1784 for hei in range(num_hei):
1785 1785 temp = jspectra[ich,:, hei]#/norm[ich,hei]
1786 1786 temp -= inttef
1787 1787 temp += jnoise[ich]
1788 1788 # print(jspectra.shape, temp.shape)
1789 1789 jspectra[ich,:, hei] = temp
1790 1790
1791 1791 # Guardar Resultados
1792 1792 self.dataOut.data_spc = jspectra
1793 1793 #self.dataOut.data_cspc = jcspectra
1794 1794
1795 1795 return 1
1796 1796
1797 1797 def removeInterference2(self):
1798 1798
1799 1799 cspc = self.dataOut.data_cspc
1800 1800 spc = self.dataOut.data_spc
1801 1801 Heights = numpy.arange(cspc.shape[2])
1802 1802 realCspc = numpy.abs(cspc)
1803 1803
1804 1804 for i in range(cspc.shape[0]):
1805 1805 LinePower= numpy.sum(realCspc[i], axis=0)
1806 1806 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1807 1807 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1808 1808 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1809 1809 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1810 1810 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1811 1811
1812 1812
1813 1813 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1814 1814 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1815 1815 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1816 1816 cspc[i,InterferenceRange,:] = numpy.NaN
1817 1817
1818 1818 self.dataOut.data_cspc = cspc
1819 1819
1820 1820 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1821 1821
1822 1822 jspectra = self.dataOut.data_spc
1823 1823 jcspectra = self.dataOut.data_cspc
1824 1824 jnoise = self.dataOut.getNoise()
1825 1825 #num_incoh = self.dataOut.nIncohInt
1826 1826 num_incoh = self.dataOut.max_nIncohInt
1827 1827 #print("spc: ", jspectra.shape, jcspectra)
1828 1828 num_channel = jspectra.shape[0]
1829 1829 num_prof = jspectra.shape[1]
1830 1830 num_hei = jspectra.shape[2]
1831 1831
1832 1832 count_hei = nhei_interf
1833 1833 # hei_interf
1834 1834 if hei_interf is None:
1835 1835 count_hei = int(num_hei / 2) # a half of total ranges
1836 1836 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1837 1837 hei_interf = numpy.asarray(hei_interf)[0]
1838 1838 #print(hei_interf)
1839 1839 # nhei_interf
1840 1840 if (nhei_interf == None):
1841 1841 nhei_interf = 5
1842 1842 if (nhei_interf < 1):
1843 1843 nhei_interf = 1
1844 1844 if (nhei_interf > count_hei):
1845 1845 nhei_interf = count_hei
1846 1846 if (offhei_interf == None):
1847 1847 offhei_interf = 0
1848 1848
1849 1849 ind_hei = list(range(num_hei))
1850 1850 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1851 1851 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1852 1852 mask_prof = numpy.asarray(list(range(num_prof)))
1853 1853 num_mask_prof = mask_prof.size
1854 1854 comp_mask_prof = [0, num_prof / 2]
1855 1855
1856 1856 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1857 1857 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1858 1858 jnoise = numpy.nan
1859 1859 noise_exist = jnoise[0] < numpy.Inf
1860 1860
1861 1861 # Subrutina de Remocion de la Interferencia
1862 1862 for ich in range(num_channel):
1863 1863 # Se ordena los espectros segun su potencia (menor a mayor)
1864 1864 power = jspectra[ich, mask_prof, :]
1865 1865 power = power[:, hei_interf]
1866 1866 power = power.sum(axis=0)
1867 1867 psort = power.ravel().argsort()
1868 1868 #print(hei_interf[psort[list(range(offhei_interf, nhei_interf + offhei_interf))]])
1869 1869 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1870 1870 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1871 1871 offhei_interf, nhei_interf + offhei_interf))]]]
1872 1872
1873 1873 if noise_exist:
1874 1874 # tmp_noise = jnoise[ich] / num_prof
1875 1875 tmp_noise = jnoise[ich]
1876 1876 junkspc_interf = junkspc_interf - tmp_noise
1877 1877 #junkspc_interf[:,comp_mask_prof] = 0
1878 1878 #print(junkspc_interf.shape)
1879 1879 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1880 1880 jspc_interf = jspc_interf.transpose()
1881 1881 # Calculando el espectro de interferencia promedio
1882 1882 noiseid = numpy.where(jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1883 1883 noiseid = noiseid[0]
1884 1884 cnoiseid = noiseid.size
1885 1885 interfid = numpy.where(jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1886 1886 interfid = interfid[0]
1887 1887 cinterfid = interfid.size
1888 1888
1889 1889 if (cnoiseid > 0):
1890 1890 jspc_interf[noiseid] = 0
1891 1891 # Expandiendo los perfiles a limpiar
1892 1892 if (cinterfid > 0):
1893 1893 new_interfid = (
1894 1894 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1895 1895 new_interfid = numpy.asarray(new_interfid)
1896 1896 new_interfid = {x for x in new_interfid}
1897 1897 new_interfid = numpy.array(list(new_interfid))
1898 1898 new_cinterfid = new_interfid.size
1899 1899 else:
1900 1900 new_cinterfid = 0
1901 1901
1902 1902 for ip in range(new_cinterfid):
1903 1903 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1904 1904 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1905 1905
1906 1906 jspectra[ich, :, ind_hei] = jspectra[ich, :,ind_hei] - jspc_interf # Corregir indices
1907 1907
1908 1908 # Removiendo la interferencia del punto de mayor interferencia
1909 1909 ListAux = jspc_interf[mask_prof].tolist()
1910 1910 maxid = ListAux.index(max(ListAux))
1911 1911 #print(cinterfid)
1912 1912 if cinterfid > 0:
1913 1913 for ip in range(cinterfid * (interf == 2) - 1):
1914 1914 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1915 1915 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1916 1916 cind = len(ind)
1917 1917
1918 1918 if (cind > 0):
1919 1919 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1920 1920 (1 + (numpy.random.uniform(cind) - 0.5) /
1921 1921 numpy.sqrt(num_incoh))
1922 1922
1923 1923 ind = numpy.array([-2, -1, 1, 2])
1924 1924 xx = numpy.zeros([4, 4])
1925 1925
1926 1926 for id1 in range(4):
1927 1927 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1928 1928 xx_inv = numpy.linalg.inv(xx)
1929 1929 xx = xx_inv[:, 0]
1930 1930 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1931 1931 yy = jspectra[ich, mask_prof[ind], :]
1932 1932 jspectra[ich, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1933 1933
1934 1934 indAux = (jspectra[ich, :, :] < tmp_noise *
1935 1935 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1936 1936 #print(indAux)
1937 1937 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1938 1938 (1 - 1 / numpy.sqrt(num_incoh))
1939 1939
1940 1940 # Remocion de Interferencia en el Cross Spectra
1941 1941 if jcspectra is None:
1942 1942 return jspectra, jcspectra
1943 1943 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1944 1944 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1945 1945
1946 1946 for ip in range(num_pairs):
1947 1947
1948 1948 #-------------------------------------------
1949 1949
1950 1950 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1951 1951 cspower = cspower[:, hei_interf]
1952 1952 cspower = cspower.sum(axis=0)
1953 1953
1954 1954 cspsort = cspower.ravel().argsort()
1955 1955 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1956 1956 offhei_interf, nhei_interf + offhei_interf))]]]
1957 1957 junkcspc_interf = junkcspc_interf.transpose()
1958 1958 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1959 1959
1960 1960 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1961 1961
1962 1962 median_real = int(numpy.median(numpy.real(
1963 1963 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1964 1964 median_imag = int(numpy.median(numpy.imag(
1965 1965 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1966 1966 comp_mask_prof = [int(e) for e in comp_mask_prof]
1967 1967 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1968 1968 median_real, median_imag)
1969 1969
1970 1970 for iprof in range(num_prof):
1971 1971 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1972 1972 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1973 1973
1974 1974 # Removiendo la Interferencia
1975 1975 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1976 1976 :, ind_hei] - jcspc_interf
1977 1977
1978 1978 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1979 1979 maxid = ListAux.index(max(ListAux))
1980 1980
1981 1981 ind = numpy.array([-2, -1, 1, 2])
1982 1982 xx = numpy.zeros([4, 4])
1983 1983
1984 1984 for id1 in range(4):
1985 1985 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1986 1986
1987 1987 xx_inv = numpy.linalg.inv(xx)
1988 1988 xx = xx_inv[:, 0]
1989 1989
1990 1990 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1991 1991 yy = jcspectra[ip, mask_prof[ind], :]
1992 1992 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1993 1993
1994 1994 # Guardar Resultados
1995 1995 self.dataOut.data_spc = jspectra
1996 1996 self.dataOut.data_cspc = jcspectra
1997 1997
1998 1998 return 1
1999 1999
2000 2000 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1, minHei=None, maxHei=None):
2001 2001
2002 2002 self.dataOut = dataOut
2003 2003
2004 2004 if mode == 1:
2005 2005 self.removeInterference(interf = 2,hei_interf = hei_interf, nhei_interf = nhei_interf, offhei_interf = offhei_interf)
2006 2006 elif mode == 2:
2007 2007 self.removeInterference2()
2008 2008 elif mode == 3:
2009 2009 self.removeInterference3(min_hei=minHei, max_hei=maxHei)
2010 2010 return self.dataOut
2011 2011
2012 2012
2013 2013 class IncohInt(Operation):
2014 2014
2015 2015 __profIndex = 0
2016 2016 __withOverapping = False
2017 2017
2018 2018 __byTime = False
2019 2019 __initime = None
2020 2020 __lastdatatime = None
2021 2021 __integrationtime = None
2022 2022
2023 2023 __buffer_spc = None
2024 2024 __buffer_cspc = None
2025 2025 __buffer_dc = None
2026 2026
2027 2027 __dataReady = False
2028 2028
2029 2029 __timeInterval = None
2030 2030 incohInt = 0
2031 2031 nOutliers = 0
2032 2032 n = None
2033 2033
2034 2034 _flagProfilesByRange = False
2035 2035 _nProfilesByRange = 0
2036 2036 def __init__(self):
2037 2037
2038 2038 Operation.__init__(self)
2039 2039
2040 2040 def setup(self, n=None, timeInterval=None, overlapping=False):
2041 2041 """
2042 2042 Set the parameters of the integration class.
2043 2043
2044 2044 Inputs:
2045 2045
2046 2046 n : Number of coherent integrations
2047 2047 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
2048 2048 overlapping :
2049 2049
2050 2050 """
2051 2051
2052 2052 self.__initime = None
2053 2053 self.__lastdatatime = 0
2054 2054
2055 2055 self.__buffer_spc = 0
2056 2056 self.__buffer_cspc = 0
2057 2057 self.__buffer_dc = 0
2058 2058
2059 2059 self.__profIndex = 0
2060 2060 self.__dataReady = False
2061 2061 self.__byTime = False
2062 2062 self.incohInt = 0
2063 2063 self.nOutliers = 0
2064 2064 if n is None and timeInterval is None:
2065 2065 raise ValueError("n or timeInterval should be specified ...")
2066 2066
2067 2067 if n is not None:
2068 2068 self.n = int(n)
2069 2069 else:
2070 2070
2071 2071 self.__integrationtime = int(timeInterval)
2072 2072 self.n = None
2073 2073 self.__byTime = True
2074 2074
2075 2075
2076 2076
2077 2077 def putData(self, data_spc, data_cspc, data_dc):
2078 2078 """
2079 2079 Add a profile to the __buffer_spc and increase in one the __profileIndex
2080 2080
2081 2081 """
2082 2082 if data_spc.all() == numpy.nan :
2083 2083 print("nan ")
2084 2084 return
2085 2085 self.__buffer_spc += data_spc
2086 2086
2087 2087 if data_cspc is None:
2088 2088 self.__buffer_cspc = None
2089 2089 else:
2090 2090 self.__buffer_cspc += data_cspc
2091 2091
2092 2092 if data_dc is None:
2093 2093 self.__buffer_dc = None
2094 2094 else:
2095 2095 self.__buffer_dc += data_dc
2096 2096
2097 2097 self.__profIndex += 1
2098 2098
2099 2099 return
2100 2100
2101 2101 def pushData(self):
2102 2102 """
2103 2103 Return the sum of the last profiles and the profiles used in the sum.
2104 2104
2105 2105 Affected:
2106 2106
2107 2107 self.__profileIndex
2108 2108
2109 2109 """
2110 2110
2111 2111 data_spc = self.__buffer_spc
2112 2112 data_cspc = self.__buffer_cspc
2113 2113 data_dc = self.__buffer_dc
2114 2114 n = self.__profIndex
2115 2115
2116 2116 self.__buffer_spc = 0
2117 2117 self.__buffer_cspc = 0
2118 2118 self.__buffer_dc = 0
2119 2119
2120 2120
2121 2121 return data_spc, data_cspc, data_dc, n
2122 2122
2123 2123 def byProfiles(self, *args):
2124 2124
2125 2125 self.__dataReady = False
2126 2126 avgdata_spc = None
2127 2127 avgdata_cspc = None
2128 2128 avgdata_dc = None
2129 2129
2130 2130 self.putData(*args)
2131 2131
2132 2132 if self.__profIndex == self.n:
2133 2133
2134 2134 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2135 2135 self.n = n
2136 2136 self.__dataReady = True
2137 2137
2138 2138 return avgdata_spc, avgdata_cspc, avgdata_dc
2139 2139
2140 2140 def byTime(self, datatime, *args):
2141 2141
2142 2142 self.__dataReady = False
2143 2143 avgdata_spc = None
2144 2144 avgdata_cspc = None
2145 2145 avgdata_dc = None
2146 2146
2147 2147 self.putData(*args)
2148 2148
2149 2149 if (datatime - self.__initime) >= self.__integrationtime:
2150 2150 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
2151 2151 self.n = n
2152 2152 self.__dataReady = True
2153 2153
2154 2154 return avgdata_spc, avgdata_cspc, avgdata_dc
2155 2155
2156 2156 def integrate(self, datatime, *args):
2157 2157
2158 2158 if self.__profIndex == 0:
2159 2159 self.__initime = datatime
2160 2160
2161 2161 if self.__byTime:
2162 2162 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
2163 2163 datatime, *args)
2164 2164 else:
2165 2165 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
2166 2166
2167 2167 if not self.__dataReady:
2168 2168 return None, None, None, None
2169 2169
2170 2170 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
2171 2171
2172 2172 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
2173 2173 if n == 1:
2174 2174 return dataOut
2175 2175
2176 2176 if dataOut.flagNoData == True:
2177 2177 return dataOut
2178 2178
2179 2179 if dataOut.flagProfilesByRange == True:
2180 2180 self._flagProfilesByRange = True
2181 2181
2182 2182 dataOut.flagNoData = True
2183 2183 dataOut.processingHeaderObj.timeIncohInt = timeInterval
2184 2184 if not self.isConfig:
2185 2185 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
2186 2186 self.setup(n, timeInterval, overlapping)
2187 2187 self.isConfig = True
2188 2188
2189 2189
2190 2190 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
2191 2191 dataOut.data_spc,
2192 2192 dataOut.data_cspc,
2193 2193 dataOut.data_dc)
2194 2194
2195 2195 self.incohInt += dataOut.nIncohInt
2196 2196
2197 2197
2198 2198 if isinstance(dataOut.data_outlier,numpy.ndarray) or isinstance(dataOut.data_outlier,int) or isinstance(dataOut.data_outlier, float):
2199 2199 self.nOutliers += dataOut.data_outlier
2200 2200
2201 2201 if self._flagProfilesByRange:
2202 2202 dataOut.flagProfilesByRange = True
2203 2203 self._nProfilesByRange += dataOut.nProfilesByRange
2204 2204
2205 2205 if self.__dataReady:
2206 print("IncohInt Done ", self.incohInt)
2207 2206 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
2208 2207 dataOut.data_spc = avgdata_spc
2209 2208 dataOut.data_cspc = avgdata_cspc
2210 2209 dataOut.data_dc = avgdata_dc
2211 2210 dataOut.nIncohInt = self.incohInt
2212 2211 dataOut.data_outlier = self.nOutliers
2213 2212 dataOut.utctime = avgdatatime
2214 2213 dataOut.flagNoData = False
2215 2214 self.incohInt = 0
2216 2215 self.nOutliers = 0
2217 2216 self.__profIndex = 0
2218 2217 dataOut.nProfilesByRange = self._nProfilesByRange
2219 2218 self._nProfilesByRange = numpy.zeros((1,len(dataOut.heightList)))
2220 2219 self._flagProfilesByRange = False
2221
2220 #print("IncohInt Done")
2222 2221 return dataOut
2223 2222
2224 2223 class dopplerFlip(Operation):
2225 2224
2226 2225 def run(self, dataOut):
2227 2226 # arreglo 1: (num_chan, num_profiles, num_heights)
2228 2227 self.dataOut = dataOut
2229 2228 # JULIA-oblicua, indice 2
2230 2229 # arreglo 2: (num_profiles, num_heights)
2231 2230 jspectra = self.dataOut.data_spc[2]
2232 2231 jspectra_tmp = numpy.zeros(jspectra.shape)
2233 2232 num_profiles = jspectra.shape[0]
2234 2233 freq_dc = int(num_profiles / 2)
2235 2234 # Flip con for
2236 2235 for j in range(num_profiles):
2237 2236 jspectra_tmp[num_profiles-j-1]= jspectra[j]
2238 2237 # Intercambio perfil de DC con perfil inmediato anterior
2239 2238 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
2240 2239 jspectra_tmp[freq_dc]= jspectra[freq_dc]
2241 2240 # canal modificado es re-escrito en el arreglo de canales
2242 2241 self.dataOut.data_spc[2] = jspectra_tmp
2243 2242
2244 2243 return self.dataOut
2245 2244
2246 2245
2247 2246
2248 2247
2249 2248
2250 2249
2251 2250 class cleanJULIAInterf(Operation):
2252 2251 """
2253 2252 OperaciΓ³n de prueba
2254 2253 """
2255 2254 __slots__ =('heights_indx', 'repeats','span' ,'step', 'factor', 'idate', 'idxs','isConfig','minHrefN', 'maxHrefN')
2256 2255 def __init__(self):
2257 2256 self.repeats = 0
2258 2257 self.factor=1
2259 2258 self.isConfig = False
2260 2259 self.idxs = None
2261 2260 self.heights_indx = None
2262 2261
2263 2262 def setup(self, dataOutHeightsList, heightsList, span=10, repeats=0, step=0, idate=None, startH=None, endH=None, minHref=None, maxHref=None):
2264 2263 totalHeihtList = dataOutHeightsList
2265 2264 heights = [float(hei) for hei in heightsList]
2266 2265 for r in range(repeats):
2267 2266 heights += [ (h+(step*(r+1))) for h in heights]
2268 2267 #print(heights)
2269 2268 self.heights_indx = [getHei_index(h,h,totalHeihtList)[0] for h in heights]
2270 2269
2271 2270 self.minHrefN, self.maxHrefN = getHei_index(minHref,maxHref,totalHeihtList)
2272 2271
2273 2272
2274 2273 self.config = True
2275 2274 self.span = span
2276 2275
2277 2276 def run(self, dataOut, heightsList, span=10, repeats=0, step=0, idate=None, startH=None, endH=None, minHref=None, maxHref=None):
2278 2277
2279 2278
2280 2279 self.dataOut = dataOut
2281 2280 startTime = datetime.datetime.combine(idate,startH)
2282 2281 endTime = datetime.datetime.combine(idate,endH)
2283 2282 currentTime = datetime.datetime.fromtimestamp(self.dataOut.utctime)
2284 2283
2285 2284 if currentTime < startTime or currentTime > endTime:
2286 2285 return self.dataOut
2287 2286
2288 2287 if not self.isConfig:
2289 2288 self.setup(self.dataOut.heightList,heightsList, span=span, repeats=repeats, step=step, idate=idate, startH=startH, endH=endH, minHref=minHref, maxHref=maxHref )
2290 2289
2291 2290 for ch in range(self.dataOut.data_spc.shape[0]):
2292 2291 i = 0
2293 2292 N_ref = self.dataOut.data_spc[ch, :, self.minHrefN: self.maxHrefN].mean()
2294 2293 mn = self.heights_indx[-1] - self.span/2
2295 2294 mx = self.heights_indx[-1] + self.span/2
2296 2295 J_lev = self.dataOut.data_spc[ch, :, mn: mx].mean() - N_ref
2297 2296
2298 2297 for hei in self.heights_indx:
2299 2298 h = hei - 1
2300 2299 mn_i = hei - self.span/2
2301 2300 mx_i = hei + self.span/2
2302 2301 self.dataOut.data_spc[ch, :,mn_i:mx_i ] -= J_lev
2303 2302 i += 1
2304 2303
2305 2304
2306 2305 return self.dataOut No newline at end of file
@@ -1,3259 +1,3249
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from schainpy.model.io.utilsIO import getHei_index
8 8 from time import time
9 9 import datetime
10 10 import numpy
11 11 #import copy
12 12 from schainpy.model.data import _noise
13 13
14 14 from matplotlib import pyplot as plt
15 15
16 16 class VoltageProc(ProcessingUnit):
17 17
18 18 def __init__(self):
19 19
20 20 ProcessingUnit.__init__(self)
21 21
22 22 self.dataOut = Voltage()
23 23 self.flip = 1
24 24 self.setupReq = False
25 25
26 26 def run(self):
27 27 #print("running volt proc")
28 28
29 29 if self.dataIn.type == 'AMISR':
30 30 self.__updateObjFromAmisrInput()
31 31
32 32 if self.dataOut.buffer_empty:
33 33 if self.dataIn.type == 'Voltage':
34 34 self.dataOut.copy(self.dataIn)
35 35 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
36 36 self.dataOut.ippSeconds = self.dataIn.ippSeconds
37 37 self.dataOut.ipp = self.dataIn.ipp
38 38
39 39 #update Processing Header:
40 40 self.dataOut.processingHeaderObj.heightList = self.dataOut.heightList
41 41 self.dataOut.processingHeaderObj.ipp = self.dataOut.ipp
42 42 self.dataOut.processingHeaderObj.nCohInt = self.dataOut.nCohInt
43 43 self.dataOut.processingHeaderObj.dtype = self.dataOut.type
44 44 self.dataOut.processingHeaderObj.channelList = self.dataOut.channelList
45 45 self.dataOut.processingHeaderObj.azimuthList = self.dataOut.azimuthList
46 46 self.dataOut.processingHeaderObj.elevationList = self.dataOut.elevationList
47 47 self.dataOut.processingHeaderObj.codeList = self.dataOut.nChannels
48 48 self.dataOut.processingHeaderObj.heightList = self.dataOut.heightList
49 49 self.dataOut.processingHeaderObj.heightResolution = self.dataOut.heightList[1] - self.dataOut.heightList[0]
50 50
51 51
52 52
53 53 def __updateObjFromAmisrInput(self):
54 54
55 55 self.dataOut.timeZone = self.dataIn.timeZone
56 56 self.dataOut.dstFlag = self.dataIn.dstFlag
57 57 self.dataOut.errorCount = self.dataIn.errorCount
58 58 self.dataOut.useLocalTime = self.dataIn.useLocalTime
59 59
60 60 self.dataOut.flagNoData = self.dataIn.flagNoData
61 61 self.dataOut.data = self.dataIn.data
62 62 self.dataOut.utctime = self.dataIn.utctime
63 63 self.dataOut.channelList = self.dataIn.channelList
64 64 #self.dataOut.timeInterval = self.dataIn.timeInterval
65 65 self.dataOut.heightList = self.dataIn.heightList
66 66 self.dataOut.nProfiles = self.dataIn.nProfiles
67 67
68 68 self.dataOut.nCohInt = self.dataIn.nCohInt
69 69 self.dataOut.ippSeconds = self.dataIn.ippSeconds
70 70 self.dataOut.frequency = self.dataIn.frequency
71 71
72 72 self.dataOut.azimuth = self.dataIn.azimuth
73 73 self.dataOut.zenith = self.dataIn.zenith
74 74
75 75 self.dataOut.beam.codeList = self.dataIn.beam.codeList
76 76 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
77 77 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
78 78
79 79
80 80 class selectChannels(Operation):
81 81
82 82 def run(self, dataOut, channelList=[]):
83 83
84 84 if isinstance(channelList, int):
85 85 channelList = [channelList]
86 86
87 87 self.channelList = channelList
88 88 if len(self.channelList) == 0:
89 89 print("Missing channelList")
90 90 return dataOut
91 91 channelIndexList = []
92 92 if not dataOut.buffer_empty: # cuando se usa proc volts como buffer de datos
93 93 return dataOut
94 94 #print("channel List: ", dataOut.channelList)
95 95 if type(dataOut.channelList) is not list: #leer array desde HDF5
96 96 try:
97 97 dataOut.channelList = dataOut.channelList.tolist()
98 98 except Exception as e:
99 99 print("Select Channels: ",e)
100 100 for channel in self.channelList:
101 101 if channel not in dataOut.channelList:
102 102 raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
103 103
104 104 index = dataOut.channelList.index(channel)
105 105 channelIndexList.append(index)
106 106
107 107 dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
108 108
109 109 #update Processing Header:
110 110 dataOut.processingHeaderObj.channelList = dataOut.channelList
111 111 dataOut.processingHeaderObj.elevationList = dataOut.elevationList
112 112 dataOut.processingHeaderObj.azimuthList = dataOut.azimuthList
113 113 dataOut.processingHeaderObj.codeList = dataOut.codeList
114 114 dataOut.processingHeaderObj.nChannels = len(dataOut.channelList)
115 115
116 116 return dataOut
117 117
118 118 def selectChannelsByIndex(self, dataOut, channelIndexList):
119 119 """
120 120 Selecciona un bloque de datos en base a canales segun el channelIndexList
121 121
122 122 Input:
123 123 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
124 124
125 125 Affected:
126 126 dataOut.data
127 127 dataOut.channelIndexList
128 128 dataOut.nChannels
129 129 dataOut.m_ProcessingHeader.totalSpectra
130 130 dataOut.systemHeaderObj.numChannels
131 131 dataOut.m_ProcessingHeader.blockSize
132 132
133 133 Return:
134 134 None
135 135 """
136 136 #print("selectChannelsByIndex")
137 137 # for channelIndex in channelIndexList:
138 138 # if channelIndex not in dataOut.channelIndexList:
139 139 # raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
140 140
141 141 if dataOut.type == 'Voltage':
142 142 if dataOut.flagDataAsBlock:
143 143 """
144 144 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
145 145 """
146 146 data = dataOut.data[channelIndexList,:,:]
147 147 else:
148 148 data = dataOut.data[channelIndexList,:]
149 149
150 150 dataOut.data = data
151 151 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
152 152 dataOut.channelList = [n for n in range(len(channelIndexList))]
153 153
154 154 elif dataOut.type == 'Spectra':
155 155 if hasattr(dataOut, 'data_spc'):
156 156 if dataOut.data_spc is None:
157 157 raise ValueError("data_spc is None")
158 158 return dataOut
159 159 else:
160 160 data_spc = dataOut.data_spc[channelIndexList, :]
161 161 dataOut.data_spc = data_spc
162 162
163 163 # if hasattr(dataOut, 'data_dc') :# and
164 164 # if dataOut.data_dc is None:
165 165 # raise ValueError("data_dc is None")
166 166 # return dataOut
167 167 # else:
168 168 # data_dc = dataOut.data_dc[channelIndexList, :]
169 169 # dataOut.data_dc = data_dc
170 170 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
171 171 dataOut.channelList = channelIndexList
172 172 dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
173 173
174 174 # channelIndexList = numpy.asarray(channelIndexList)
175 175 dataOut.elevationList = numpy.asarray(dataOut.elevationList)
176 176 dataOut.azimuthList = numpy.asarray(dataOut.azimuthList)
177 177 dataOut.codeList = numpy.asarray(dataOut.codeList)
178 178 if (len(dataOut.elevationList) > 0):
179 179 dataOut.elevationList = dataOut.elevationList[channelIndexList]
180 180 dataOut.azimuthList = dataOut.azimuthList[channelIndexList]
181 181 dataOut.codeList = dataOut.codeList[channelIndexList]
182 182
183 183 return dataOut
184 184
185 185 def __selectPairsByChannel(self, dataOut, channelList=None):
186 186 #print("__selectPairsByChannel")
187 187 if channelList == None:
188 188 return
189 189
190 190 pairsIndexListSelected = []
191 191 for pairIndex in dataOut.pairsIndexList:
192 192 # First pair
193 193 if dataOut.pairsList[pairIndex][0] not in channelList:
194 194 continue
195 195 # Second pair
196 196 if dataOut.pairsList[pairIndex][1] not in channelList:
197 197 continue
198 198
199 199 pairsIndexListSelected.append(pairIndex)
200 200 if not pairsIndexListSelected:
201 201 dataOut.data_cspc = None
202 202 dataOut.pairsList = []
203 203 return
204 204
205 205 dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
206 206 dataOut.pairsList = [dataOut.pairsList[i]
207 207 for i in pairsIndexListSelected]
208 208
209 209 return dataOut
210 210
211 211 class selectHeights(Operation):
212 212
213 213 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
214 214 """
215 215 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
216 216 minHei <= height <= maxHei
217 217
218 218 Input:
219 219 minHei : valor minimo de altura a considerar
220 220 maxHei : valor maximo de altura a considerar
221 221
222 222 Affected:
223 223 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
224 224
225 225 Return:
226 226 1 si el metodo se ejecuto con exito caso contrario devuelve 0
227 227 """
228 228
229 229 self.dataOut = dataOut
230 230
231 231 if minHei and maxHei:
232 232
233 233 if (minHei < dataOut.heightList[0]):
234 234 minHei = dataOut.heightList[0]
235 235
236 236 if (maxHei > dataOut.heightList[-1]):
237 237 maxHei = dataOut.heightList[-1]
238 238
239 239 minIndex = 0
240 240 maxIndex = 0
241 241 heights = dataOut.heightList
242 242
243 243 inda = numpy.where(heights >= minHei)
244 244 indb = numpy.where(heights <= maxHei)
245 245
246 246 try:
247 247 minIndex = inda[0][0]
248 248 except:
249 249 minIndex = 0
250 250
251 251 try:
252 252 maxIndex = indb[0][-1]
253 253 except:
254 254 maxIndex = len(heights)
255 255
256 256 self.selectHeightsByIndex(minIndex, maxIndex)
257 257
258 258 #update Processing Header:
259 259 dataOut.processingHeaderObj.heightList = dataOut.heightList
260 260
261 261
262 262
263 263 return dataOut
264 264
265 265 def selectHeightsByIndex(self, minIndex, maxIndex):
266 266 """
267 267 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
268 268 minIndex <= index <= maxIndex
269 269
270 270 Input:
271 271 minIndex : valor de indice minimo de altura a considerar
272 272 maxIndex : valor de indice maximo de altura a considerar
273 273
274 274 Affected:
275 275 self.dataOut.data
276 276 self.dataOut.heightList
277 277
278 278 Return:
279 279 1 si el metodo se ejecuto con exito caso contrario devuelve 0
280 280 """
281 281
282 282 if self.dataOut.type == 'Voltage':
283 283 if (minIndex < 0) or (minIndex > maxIndex):
284 284 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
285 285
286 286 if (maxIndex >= self.dataOut.nHeights):
287 287 maxIndex = self.dataOut.nHeights
288 288
289 289 #voltage
290 290 if self.dataOut.flagDataAsBlock:
291 291 """
292 292 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
293 293 """
294 294 data = self.dataOut.data[:,:, minIndex:maxIndex]
295 295 else:
296 296 data = self.dataOut.data[:, minIndex:maxIndex]
297 297
298 298 # firstHeight = self.dataOut.heightList[minIndex]
299 299
300 300 self.dataOut.data = data
301 301 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
302 302
303 303 if self.dataOut.nHeights <= 1:
304 304 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
305 305 elif self.dataOut.type == 'Spectra':
306 306 if (minIndex < 0) or (minIndex > maxIndex):
307 307 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
308 308 minIndex, maxIndex))
309 309
310 310 if (maxIndex >= self.dataOut.nHeights):
311 311 maxIndex = self.dataOut.nHeights - 1
312 312
313 313 # Spectra
314 314 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
315 315
316 316 data_cspc = None
317 317 if self.dataOut.data_cspc is not None:
318 318 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
319 319
320 320 data_dc = None
321 321 if self.dataOut.data_dc is not None:
322 322 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
323 323
324 324 self.dataOut.data_spc = data_spc
325 325 self.dataOut.data_cspc = data_cspc
326 326 self.dataOut.data_dc = data_dc
327 327
328 328 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
329 329
330 330 return 1
331 331
332 332
333 333 class filterByHeights(Operation):
334 334 ifConfig=False
335 335 deltaHeight = None
336 336 newdelta=None
337 337 newheights=None
338 338 r=None
339 339 h0=None
340 340 nHeights=None
341 341 def run(self, dataOut, window):
342 342
343 343
344 344 # print("1",dataOut.data.shape)
345 345 # print(dataOut.nHeights)
346 346 if window == None:
347 347 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / self.deltaHeight
348 348
349 349 if not self.ifConfig: #and dataOut.useInputBuffer:
350 350 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
351 351 self.ifConfig = True
352 352 self.newdelta = self.deltaHeight * window
353 353 self.r = dataOut.nHeights % window
354 354 self.newheights = (dataOut.nHeights-self.r)/window
355 355 self.h0 = dataOut.heightList[0]
356 356 self.nHeights = dataOut.nHeights
357 357 if self.newheights <= 1:
358 358 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
359 359
360 360 if dataOut.flagDataAsBlock:
361 361 """
362 362 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
363 363 """
364 364 buffer = dataOut.data[:, :, 0:int(self.nHeights-self.r)]
365 365 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(self.nHeights/window), window)
366 366 buffer = numpy.sum(buffer,3)
367 367
368 368 else:
369 369 buffer = dataOut.data[:,0:int(self.nHeights-self.r)]
370 370 buffer = buffer.reshape(dataOut.nChannels,int(self.nHeights/window),int(window))
371 371 buffer = numpy.sum(buffer,2)
372 372
373 373 dataOut.data = buffer
374 374 dataOut.heightList = self.h0 + numpy.arange( self.newheights )*self.newdelta
375 375 dataOut.windowOfFilter = window
376 376
377 377 #update Processing Header:
378 378 dataOut.processingHeaderObj.heightList = dataOut.heightList
379 379 dataOut.processingHeaderObj.nWindows = window
380 380
381 381 return dataOut
382 382
383 383
384 384
385 385 class setH0(Operation):
386 386
387 387 def run(self, dataOut, h0, deltaHeight = None):
388 388
389 389 if not deltaHeight:
390 390 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
391 391
392 392 nHeights = dataOut.nHeights
393 393
394 394 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
395 395
396 396 dataOut.heightList = newHeiRange
397 397
398 398 #update Processing Header:
399 399 dataOut.processingHeaderObj.heightList = dataOut.heightList
400 400
401 401 return dataOut
402 402
403 403
404 404 class deFlip(Operation):
405 405
406 406 def run(self, dataOut, channelList = []):
407 407
408 408 data = dataOut.data.copy()
409 409
410 410 if dataOut.flagDataAsBlock:
411 411 flip = self.flip
412 412 profileList = list(range(dataOut.nProfiles))
413 413
414 414 if not channelList:
415 415 for thisProfile in profileList:
416 416 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
417 417 flip *= -1.0
418 418 else:
419 419 for thisChannel in channelList:
420 420 if thisChannel not in dataOut.channelList:
421 421 continue
422 422
423 423 for thisProfile in profileList:
424 424 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
425 425 flip *= -1.0
426 426
427 427 self.flip = flip
428 428
429 429 else:
430 430 if not channelList:
431 431 data[:,:] = data[:,:]*self.flip
432 432 else:
433 433 for thisChannel in channelList:
434 434 if thisChannel not in dataOut.channelList:
435 435 continue
436 436
437 437 data[thisChannel,:] = data[thisChannel,:]*self.flip
438 438
439 439 self.flip *= -1.
440 440
441 441 dataOut.data = data
442 442
443 443 return dataOut
444 444
445 445
446 446 class setAttribute(Operation):
447 447 '''
448 448 Set an arbitrary attribute(s) to dataOut
449 449 '''
450 450
451 451 def __init__(self):
452 452
453 453 Operation.__init__(self)
454 454 self._ready = False
455 455
456 456 def run(self, dataOut, **kwargs):
457 457
458 458 for key, value in kwargs.items():
459 459 setattr(dataOut, key, value)
460 460
461 461 return dataOut
462 462
463 463
464 464 @MPDecorator
465 465 class printAttribute(Operation):
466 466 '''
467 467 Print an arbitrary attribute of dataOut
468 468 '''
469 469
470 470 def __init__(self):
471 471
472 472 Operation.__init__(self)
473 473
474 474 def run(self, dataOut, attributes):
475 475
476 476 if isinstance(attributes, str):
477 477 attributes = [attributes]
478 478 for attr in attributes:
479 479 if hasattr(dataOut, attr):
480 480 log.log(getattr(dataOut, attr), attr)
481 481
482 482 class cleanHeightsInterf(Operation):
483 483 __slots__ =('heights_indx', 'repeats', 'step', 'factor', 'idate', 'idxs','config','wMask')
484 484 def __init__(self):
485 485 self.repeats = 0
486 486 self.factor=1
487 487 self.wMask = None
488 488 self.config = False
489 489 self.idxs = None
490 490 self.heights_indx = None
491 491
492 492 def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None):
493 493
494 494 #print(dataOut.data.shape)
495 495
496 496 startTime = datetime.datetime.combine(idate,startH)
497 497 endTime = datetime.datetime.combine(idate,endH)
498 498 currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
499 499
500 500 if currentTime < startTime or currentTime > endTime:
501 501 return dataOut
502 502 if not self.config:
503 503
504 504 #print(wMask)
505 505 heights = [float(hei) for hei in heightsList]
506 506 for r in range(repeats):
507 507 heights += [ (h+(step*(r+1))) for h in heights]
508 508 #print(heights)
509 509 heiList = dataOut.heightList
510 510 self.heights_indx = [getHei_index(h,h,heiList)[0] for h in heights]
511 511
512 512 self.wMask = numpy.asarray(factor)
513 513 self.wMask = numpy.tile(self.wMask,(repeats+2))
514 514 self.config = True
515 515
516 516 """
517 517 getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None)
518 518 """
519 519 #print(self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref)))
520 520
521 521
522 522 for ch in range(dataOut.data.shape[0]):
523 523 i = 0
524 524
525 525
526 526 for hei in self.heights_indx:
527 527 h = hei - 1
528 528
529 529
530 530 if dataOut.data.ndim < 3:
531 531 module = numpy.absolute(dataOut.data[ch,h])
532 532 prev_h1 = numpy.absolute(dataOut.data[ch,h-1])
533 533 dataOut.data[ch,h] = (dataOut.data[ch,h])/module * prev_h1
534 534
535 535 #dataOut.data[ch,hei-1] = (dataOut.data[ch,hei-1])*self.wMask[i]
536 536 else:
537 537 module = numpy.absolute(dataOut.data[ch,:,h])
538 538 prev_h1 = numpy.absolute(dataOut.data[ch,:,h-1])
539 539 dataOut.data[ch,:,h] = (dataOut.data[ch,:,h])/module * prev_h1
540 540 #dataOut.data[ch,:,hei-1] = (dataOut.data[ch,:,hei-1])*self.wMask[i]
541 541 #print("done")
542 542 i += 1
543 543
544 544
545 545 return dataOut
546 546
547 547
548 548
549 549 class interpolateHeights(Operation):
550 550
551 551 def run(self, dataOut, topLim, botLim):
552 552 #69 al 72 para julia
553 553 #82-84 para meteoros
554 554 if len(numpy.shape(dataOut.data))==2:
555 555 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
556 556 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
557 557 #dataOut.data[:,botLim:limSup+1] = sampInterp
558 558 dataOut.data[:,botLim:topLim+1] = sampInterp
559 559 else:
560 560 nHeights = dataOut.data.shape[2]
561 561 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
562 562 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
563 563 f = interpolate.interp1d(x, y, axis = 2)
564 564 xnew = numpy.arange(botLim,topLim+1)
565 565 ynew = f(xnew)
566 566 dataOut.data[:,:,botLim:topLim+1] = ynew
567 567
568 568 return dataOut
569 569
570 570
571 571 class CohInt(Operation):
572 572
573 573 isConfig = False
574 574 __profIndex = 0
575 575 __byTime = False
576 576 __initime = None
577 577 __lastdatatime = None
578 578 __integrationtime = None
579 579 __buffer = None
580 580 __bufferStride = []
581 581 __dataReady = False
582 582 __profIndexStride = 0
583 583 __dataToPutStride = False
584 584 n = None
585 585
586 586 def __init__(self, **kwargs):
587 587
588 588 Operation.__init__(self, **kwargs)
589 589
590 590 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
591 591 """
592 592 Set the parameters of the integration class.
593 593
594 594 Inputs:
595 595
596 596 n : Number of coherent integrations
597 597 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
598 598 overlapping :
599 599 """
600 600
601 601 self.__initime = None
602 602 self.__lastdatatime = 0
603 603 self.__buffer = None
604 604 self.__dataReady = False
605 605 self.byblock = byblock
606 606 self.stride = stride
607 607
608 608 if n == None and timeInterval == None:
609 609 raise ValueError("n or timeInterval should be specified ...")
610 610
611 611 if n != None:
612 612 self.n = n
613 613 self.__byTime = False
614 614 else:
615 615 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
616 616 self.n = 9999
617 617 self.__byTime = True
618 618
619 619 if overlapping:
620 620 self.__withOverlapping = True
621 621 self.__buffer = None
622 622 else:
623 623 self.__withOverlapping = False
624 624 self.__buffer = 0
625 625
626 626 self.__profIndex = 0
627 627
628 628 def putData(self, data):
629 629
630 630 """
631 631 Add a profile to the __buffer and increase in one the __profileIndex
632 632
633 633 """
634 634
635 635 if not self.__withOverlapping:
636 636 self.__buffer += data.copy()
637 637 self.__profIndex += 1
638 638 return
639 639
640 640 #Overlapping data
641 641 nChannels, nHeis = data.shape
642 642 data = numpy.reshape(data, (1, nChannels, nHeis))
643 643
644 644 #If the buffer is empty then it takes the data value
645 645 if self.__buffer is None:
646 646 self.__buffer = data
647 647 self.__profIndex += 1
648 648 return
649 649
650 650 #If the buffer length is lower than n then stakcing the data value
651 651 if self.__profIndex < self.n:
652 652 self.__buffer = numpy.vstack((self.__buffer, data))
653 653 self.__profIndex += 1
654 654 return
655 655
656 656 #If the buffer length is equal to n then replacing the last buffer value with the data value
657 657 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
658 658 self.__buffer[self.n-1] = data
659 659 self.__profIndex = self.n
660 660 return
661 661
662 662
663 663 def pushData(self):
664 664 """
665 665 Return the sum of the last profiles and the profiles used in the sum.
666 666
667 667 Affected:
668 668
669 669 self.__profileIndex
670 670
671 671 """
672 672
673 673 if not self.__withOverlapping:
674 674 data = self.__buffer
675 675 n = self.__profIndex
676 676
677 677 self.__buffer = 0
678 678 self.__profIndex = 0
679 679
680 680 return data, n
681 681
682 682 #Integration with Overlapping
683 683 data = numpy.sum(self.__buffer, axis=0)
684 684 # print data
685 685 # raise
686 686 n = self.__profIndex
687 687
688 688 return data, n
689 689
690 690 def byProfiles(self, data):
691 691
692 692 self.__dataReady = False
693 693 avgdata = None
694 694 # n = None
695 695 # print data
696 696 # raise
697 697 self.putData(data)
698 698
699 699 if self.__profIndex == self.n:
700 700 avgdata, n = self.pushData()
701 701 self.__dataReady = True
702 702
703 703 return avgdata
704 704
705 705 def byTime(self, data, datatime):
706 706
707 707 self.__dataReady = False
708 708 avgdata = None
709 709 n = None
710 710
711 711 self.putData(data)
712 712
713 713 if (datatime - self.__initime) >= self.__integrationtime:
714 714 avgdata, n = self.pushData()
715 715 self.n = n
716 716 self.__dataReady = True
717 717
718 718 return avgdata
719 719
720 720 def integrateByStride(self, data, datatime):
721 721 # print data
722 722 if self.__profIndex == 0:
723 723 self.__buffer = [[data.copy(), datatime]]
724 724 else:
725 725 self.__buffer.append([data.copy(),datatime])
726 726 self.__profIndex += 1
727 727 self.__dataReady = False
728 728
729 729 if self.__profIndex == self.n * self.stride :
730 730 self.__dataToPutStride = True
731 731 self.__profIndexStride = 0
732 732 self.__profIndex = 0
733 733 self.__bufferStride = []
734 734 for i in range(self.stride):
735 735 current = self.__buffer[i::self.stride]
736 736 data = numpy.sum([t[0] for t in current], axis=0)
737 737 avgdatatime = numpy.average([t[1] for t in current])
738 738 # print data
739 739 self.__bufferStride.append((data, avgdatatime))
740 740
741 741 if self.__dataToPutStride:
742 742 self.__dataReady = True
743 743 self.__profIndexStride += 1
744 744 if self.__profIndexStride == self.stride:
745 745 self.__dataToPutStride = False
746 746 # print self.__bufferStride[self.__profIndexStride - 1]
747 747 # raise
748 748 return self.__bufferStride[self.__profIndexStride - 1]
749 749
750 750
751 751 return None, None
752 752
753 753 def integrate(self, data, datatime=None):
754 754
755 755 if self.__initime == None:
756 756 self.__initime = datatime
757 757
758 758 if self.__byTime:
759 759 avgdata = self.byTime(data, datatime)
760 760 else:
761 761 avgdata = self.byProfiles(data)
762 762
763 763
764 764 self.__lastdatatime = datatime
765 765
766 766 if avgdata is None:
767 767 return None, None
768 768
769 769 avgdatatime = self.__initime
770 770
771 771 deltatime = datatime - self.__lastdatatime
772 772
773 773 if not self.__withOverlapping:
774 774 self.__initime = datatime
775 775 else:
776 776 self.__initime += deltatime
777 777
778 778 return avgdata, avgdatatime
779 779
780 780 def integrateByBlock(self, dataOut):
781 781
782 782 times = int(dataOut.data.shape[1]/self.n)
783 783 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
784 784
785 785 id_min = 0
786 786 id_max = self.n
787 787
788 788 for i in range(times):
789 789 junk = dataOut.data[:,id_min:id_max,:]
790 790 avgdata[:,i,:] = junk.sum(axis=1)
791 791 id_min += self.n
792 792 id_max += self.n
793 793
794 794 timeInterval = dataOut.ippSeconds*self.n
795 795 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
796 796 self.__dataReady = True
797 797 return avgdata, avgdatatime
798 798
799 799 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
800 800
801 801 if not self.isConfig:
802 802 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
803 803 self.isConfig = True
804 804
805 805 if dataOut.flagDataAsBlock:
806 806 """
807 807 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
808 808 """
809 809 avgdata, avgdatatime = self.integrateByBlock(dataOut)
810 810 dataOut.nProfiles /= self.n
811 811 else:
812 812 if stride is None:
813 813 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
814 814 else:
815 815 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
816 816
817 817
818 818 # dataOut.timeInterval *= n
819 819 dataOut.flagNoData = True
820 820
821 821 if self.__dataReady:
822 822 dataOut.data = avgdata
823 823 if not dataOut.flagCohInt:
824 824 dataOut.nCohInt *= self.n
825 825 dataOut.flagCohInt = True
826 826 dataOut.utctime = avgdatatime
827 827 # print avgdata, avgdatatime
828 828 # raise
829 829 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
830 830 dataOut.flagNoData = False
831 831
832 832 #update Processing Header:
833 833 dataOut.processingHeaderObj.nCohInt = dataOut.nCohInt
834 834
835 835
836 836 return dataOut
837 837
838 838 class Decoder(Operation):
839 839
840 840 isConfig = False
841 841 __profIndex = 0
842 842
843 843 code = None
844 844
845 845 nCode = None
846 846 nBaud = None
847 847
848 848 def __init__(self, **kwargs):
849 849
850 850 Operation.__init__(self, **kwargs)
851 851
852 852 self.times = None
853 853 self.osamp = None
854 854 # self.__setValues = False
855 855 self.isConfig = False
856 856 self.setupReq = False
857 857 def setup(self, code, osamp, dataOut):
858 858
859 859 self.__profIndex = 0
860 860
861 861 self.code = code
862 862
863 863 self.nCode = len(code)
864 864 self.nBaud = len(code[0])
865 865 if (osamp != None) and (osamp >1):
866 866 self.osamp = osamp
867 867 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
868 868 self.nBaud = self.nBaud*self.osamp
869 869
870 870 self.__nChannels = dataOut.nChannels
871 871 self.__nProfiles = dataOut.nProfiles
872 872 self.__nHeis = dataOut.nHeights
873 873
874 874 if self.__nHeis < self.nBaud:
875 875 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
876 876
877 877 #Frequency
878 878 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
879 879
880 880 __codeBuffer[:,0:self.nBaud] = self.code
881 881
882 882 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
883 883
884 884 if dataOut.flagDataAsBlock:
885 885
886 886 self.ndatadec = self.__nHeis #- self.nBaud + 1
887 887
888 888 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
889 889
890 890 else:
891 891
892 892 #Time
893 893 self.ndatadec = self.__nHeis #- self.nBaud + 1
894 894
895 895 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
896 896
897 897 def __convolutionInFreq(self, data):
898 898
899 899 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
900 900
901 901 fft_data = numpy.fft.fft(data, axis=1)
902 902
903 903 conv = fft_data*fft_code
904 904
905 905 data = numpy.fft.ifft(conv,axis=1)
906 906
907 907 return data
908 908
909 909 def __convolutionInFreqOpt(self, data):
910 910
911 911 raise NotImplementedError
912 912
913 913 def __convolutionInTime(self, data):
914 914
915 915 code = self.code[self.__profIndex]
916 916 for i in range(self.__nChannels):
917 917 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
918 918
919 919 return self.datadecTime
920 920
921 921 def __convolutionByBlockInTime(self, data):
922 922
923 923 repetitions = int(self.__nProfiles / self.nCode)
924 924 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
925 925 junk = junk.flatten()
926 926 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
927 927 profilesList = range(self.__nProfiles)
928 928
929 929 for i in range(self.__nChannels):
930 930 for j in profilesList:
931 931 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
932 932 return self.datadecTime
933 933
934 934 def __convolutionByBlockInFreq(self, data):
935 935
936 936 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
937 937
938 938
939 939 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
940 940
941 941 fft_data = numpy.fft.fft(data, axis=2)
942 942
943 943 conv = fft_data*fft_code
944 944
945 945 data = numpy.fft.ifft(conv,axis=2)
946 946
947 947 return data
948 948
949 949
950 950 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
951 951
952 952 if dataOut.flagDecodeData:
953 953 print("This data is already decoded, recoding again ...")
954 954
955 955 if not self.isConfig:
956 956
957 957 if code is None:
958 958 if dataOut.code is None:
959 959 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
960 960
961 961 code = dataOut.code
962 962 else:
963 963 code = numpy.array(code).reshape(nCode,nBaud)
964 964 self.setup(code, osamp, dataOut)
965 965
966 966 self.isConfig = True
967 967
968 968 if mode == 3:
969 969 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
970 970
971 971 if times != None:
972 972 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
973 973
974 974 if self.code is None:
975 975 print("Fail decoding: Code is not defined.")
976 976 return
977 977
978 978 self.__nProfiles = dataOut.nProfiles
979 979 datadec = None
980 980
981 981 if mode == 3:
982 982 mode = 0
983 983
984 984 if dataOut.flagDataAsBlock:
985 985 """
986 986 Decoding when data have been read as block,
987 987 """
988 988
989 989 if mode == 0:
990 990 datadec = self.__convolutionByBlockInTime(dataOut.data)
991 991 if mode == 1:
992 992 datadec = self.__convolutionByBlockInFreq(dataOut.data)
993 993 else:
994 994 """
995 995 Decoding when data have been read profile by profile
996 996 """
997 997 if mode == 0:
998 998 datadec = self.__convolutionInTime(dataOut.data)
999 999
1000 1000 if mode == 1:
1001 1001 datadec = self.__convolutionInFreq(dataOut.data)
1002 1002
1003 1003 if mode == 2:
1004 1004 datadec = self.__convolutionInFreqOpt(dataOut.data)
1005 1005
1006 1006 if datadec is None:
1007 1007 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
1008 1008
1009 1009 dataOut.code = self.code
1010 1010 dataOut.nCode = self.nCode
1011 1011 dataOut.nBaud = self.nBaud
1012 1012
1013 1013 dataOut.data = datadec
1014 1014 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
1015 1015 dataOut.flagDecodeData = True #asumo q la data esta decodificada
1016 1016
1017 1017
1018 1018 #update Processing Header:
1019 1019 dataOut.radarControllerHeaderObj.code = self.code
1020 1020 dataOut.radarControllerHeaderObj.nCode = self.nCode
1021 1021 dataOut.radarControllerHeaderObj.nBaud = self.nBaud
1022 1022 dataOut.radarControllerHeaderObj.nOsamp = osamp
1023 1023 #update Processing Header:
1024 1024 dataOut.processingHeaderObj.heightList = dataOut.heightList
1025 1025 dataOut.processingHeaderObj.heightResolution = dataOut.heightList[1]-dataOut.heightList[0]
1026 1026
1027 1027 if self.__profIndex == self.nCode-1:
1028 1028 self.__profIndex = 0
1029 1029 return dataOut
1030 1030
1031 1031 self.__profIndex += 1
1032 1032
1033 1033 return dataOut
1034 1034
1035 1035 class ProfileConcat(Operation):
1036 1036
1037 1037 isConfig = False
1038 1038 buffer = None
1039 1039
1040 1040 def __init__(self, **kwargs):
1041 1041
1042 1042 Operation.__init__(self, **kwargs)
1043 1043 self.profileIndex = 0
1044 1044
1045 1045 def reset(self):
1046 1046 self.buffer = numpy.zeros_like(self.buffer)
1047 1047 self.start_index = 0
1048 1048 self.times = 1
1049 1049
1050 1050 def setup(self, data, m, n=1):
1051 1051 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
1052 1052 self.nHeights = data.shape[1]#.nHeights
1053 1053 self.start_index = 0
1054 1054 self.times = 1
1055 1055
1056 1056 def concat(self, data):
1057 1057
1058 1058 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
1059 1059 self.start_index = self.start_index + self.nHeights
1060 1060
1061 1061 def run(self, dataOut, m):
1062 1062 dataOut.flagNoData = True
1063 1063
1064 1064 if not self.isConfig:
1065 1065 self.setup(dataOut.data, m, 1)
1066 1066 self.isConfig = True
1067 1067
1068 1068 if dataOut.flagDataAsBlock:
1069 1069 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
1070 1070
1071 1071 else:
1072 1072 self.concat(dataOut.data)
1073 1073 self.times += 1
1074 1074 if self.times > m:
1075 1075 dataOut.data = self.buffer
1076 1076 self.reset()
1077 1077 dataOut.flagNoData = False
1078 1078 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
1079 1079 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1080 1080 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
1081 1081 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
1082 1082 dataOut.ippSeconds *= m
1083 1083
1084 1084 #update Processing Header:
1085 1085 dataOut.processingHeaderObj.heightList = dataOut.heightList
1086 1086 dataOut.processingHeaderObj.ipp = dataOut.ippSeconds
1087 1087
1088 1088 return dataOut
1089 1089
1090 1090 class ProfileSelector(Operation):
1091 1091
1092 1092 profileIndex = None
1093 1093 # Tamanho total de los perfiles
1094 1094 nProfiles = None
1095 1095
1096 1096 def __init__(self, **kwargs):
1097 1097
1098 1098 Operation.__init__(self, **kwargs)
1099 1099 self.profileIndex = 0
1100 1100
1101 1101 def incProfileIndex(self):
1102 1102
1103 1103 self.profileIndex += 1
1104 1104
1105 1105 if self.profileIndex >= self.nProfiles:
1106 1106 self.profileIndex = 0
1107 1107
1108 1108 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
1109 1109
1110 1110 if profileIndex < minIndex:
1111 1111 return False
1112 1112
1113 1113 if profileIndex > maxIndex:
1114 1114 return False
1115 1115
1116 1116 return True
1117 1117
1118 1118 def isThisProfileInList(self, profileIndex, profileList):
1119 1119
1120 1120 if profileIndex not in profileList:
1121 1121 return False
1122 1122
1123 1123 return True
1124 1124
1125 1125 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
1126 1126
1127 1127 """
1128 1128 ProfileSelector:
1129 1129
1130 1130 Inputs:
1131 1131 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
1132 1132
1133 1133 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
1134 1134
1135 1135 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
1136 1136
1137 1137 """
1138 1138
1139 1139 if rangeList is not None:
1140 1140 if type(rangeList[0]) not in (tuple, list):
1141 1141 rangeList = [rangeList]
1142 1142
1143 1143 dataOut.flagNoData = True
1144 1144
1145 1145 if dataOut.flagDataAsBlock:
1146 1146 """
1147 1147 data dimension = [nChannels, nProfiles, nHeis]
1148 1148 """
1149 1149 if profileList != None:
1150 1150 dataOut.data = dataOut.data[:,profileList,:]
1151 1151
1152 1152 if profileRangeList != None:
1153 1153 minIndex = profileRangeList[0]
1154 1154 maxIndex = profileRangeList[1]
1155 1155 profileList = list(range(minIndex, maxIndex+1))
1156 1156
1157 1157 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
1158 1158
1159 1159 if rangeList != None:
1160 1160
1161 1161 profileList = []
1162 1162
1163 1163 for thisRange in rangeList:
1164 1164 minIndex = thisRange[0]
1165 1165 maxIndex = thisRange[1]
1166 1166
1167 1167 profileList.extend(list(range(minIndex, maxIndex+1)))
1168 1168
1169 1169 dataOut.data = dataOut.data[:,profileList,:]
1170 1170
1171 1171 dataOut.nProfiles = len(profileList)
1172 1172 dataOut.profileIndex = dataOut.nProfiles - 1
1173 1173 dataOut.flagNoData = False
1174 1174
1175 1175 return dataOut
1176 1176
1177 1177 """
1178 1178 data dimension = [nChannels, nHeis]
1179 1179 """
1180 1180
1181 1181 if profileList != None:
1182 1182
1183 1183 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1184 1184
1185 1185 self.nProfiles = len(profileList)
1186 1186 dataOut.nProfiles = self.nProfiles
1187 1187 dataOut.profileIndex = self.profileIndex
1188 1188 dataOut.flagNoData = False
1189 1189
1190 1190 self.incProfileIndex()
1191 1191 return dataOut
1192 1192
1193 1193 if profileRangeList != None:
1194 1194
1195 1195 minIndex = profileRangeList[0]
1196 1196 maxIndex = profileRangeList[1]
1197 1197
1198 1198 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1199 1199
1200 1200 self.nProfiles = maxIndex - minIndex + 1
1201 1201 dataOut.nProfiles = self.nProfiles
1202 1202 dataOut.profileIndex = self.profileIndex
1203 1203 dataOut.flagNoData = False
1204 1204
1205 1205 self.incProfileIndex()
1206 1206 return dataOut
1207 1207
1208 1208 if rangeList != None:
1209 1209
1210 1210 nProfiles = 0
1211 1211
1212 1212 for thisRange in rangeList:
1213 1213 minIndex = thisRange[0]
1214 1214 maxIndex = thisRange[1]
1215 1215
1216 1216 nProfiles += maxIndex - minIndex + 1
1217 1217
1218 1218 for thisRange in rangeList:
1219 1219
1220 1220 minIndex = thisRange[0]
1221 1221 maxIndex = thisRange[1]
1222 1222
1223 1223 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1224 1224
1225 1225 self.nProfiles = nProfiles
1226 1226 dataOut.nProfiles = self.nProfiles
1227 1227 dataOut.profileIndex = self.profileIndex
1228 1228 dataOut.flagNoData = False
1229 1229
1230 1230 self.incProfileIndex()
1231 1231
1232 1232 break
1233 1233
1234 1234 return dataOut
1235 1235
1236 1236
1237 1237 if beam != None: #beam is only for AMISR data
1238 1238 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1239 1239 dataOut.flagNoData = False
1240 1240 dataOut.profileIndex = self.profileIndex
1241 1241
1242 1242 self.incProfileIndex()
1243 1243
1244 1244 return dataOut
1245 1245
1246 1246 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1247 1247
1248 1248
1249 1249 class Reshaper(Operation):
1250 1250
1251 1251 def __init__(self, **kwargs):
1252 1252
1253 1253 Operation.__init__(self, **kwargs)
1254 1254
1255 1255 self.__buffer = None
1256 1256 self.__nitems = 0
1257 1257
1258 1258 def __appendProfile(self, dataOut, nTxs):
1259 1259
1260 1260 if self.__buffer is None:
1261 1261 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1262 1262 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1263 1263
1264 1264 ini = dataOut.nHeights * self.__nitems
1265 1265 end = ini + dataOut.nHeights
1266 1266
1267 1267 self.__buffer[:, ini:end] = dataOut.data
1268 1268
1269 1269 self.__nitems += 1
1270 1270
1271 1271 return int(self.__nitems*nTxs)
1272 1272
1273 1273 def __getBuffer(self):
1274 1274
1275 1275 if self.__nitems == int(1./self.__nTxs):
1276 1276
1277 1277 self.__nitems = 0
1278 1278
1279 1279 return self.__buffer.copy()
1280 1280
1281 1281 return None
1282 1282
1283 1283 def __checkInputs(self, dataOut, shape, nTxs):
1284 1284
1285 1285 if shape is None and nTxs is None:
1286 1286 raise ValueError("Reshaper: shape of factor should be defined")
1287 1287
1288 1288 if nTxs:
1289 1289 if nTxs < 0:
1290 1290 raise ValueError("nTxs should be greater than 0")
1291 1291
1292 1292 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1293 1293 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1294 1294
1295 1295 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1296 1296
1297 1297 return shape, nTxs
1298 1298
1299 1299 if len(shape) != 2 and len(shape) != 3:
1300 1300 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1301 1301
1302 1302 if len(shape) == 2:
1303 1303 shape_tuple = [dataOut.nChannels]
1304 1304 shape_tuple.extend(shape)
1305 1305 else:
1306 1306 shape_tuple = list(shape)
1307 1307
1308 1308 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1309 1309
1310 1310 return shape_tuple, nTxs
1311 1311
1312 1312 def run(self, dataOut, shape=None, nTxs=None):
1313 1313
1314 1314 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1315 1315
1316 1316 dataOut.flagNoData = True
1317 1317 profileIndex = None
1318 1318
1319 1319 if dataOut.flagDataAsBlock:
1320 1320
1321 1321 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1322 1322 dataOut.flagNoData = False
1323 1323
1324 1324 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1325 1325
1326 1326 else:
1327 1327
1328 1328 if self.__nTxs < 1:
1329 1329
1330 1330 self.__appendProfile(dataOut, self.__nTxs)
1331 1331 new_data = self.__getBuffer()
1332 1332
1333 1333 if new_data is not None:
1334 1334 dataOut.data = new_data
1335 1335 dataOut.flagNoData = False
1336 1336
1337 1337 profileIndex = dataOut.profileIndex*nTxs
1338 1338
1339 1339 else:
1340 1340 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1341 1341
1342 1342 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1343 1343
1344 1344 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1345 1345
1346 1346 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1347 1347
1348 1348 dataOut.profileIndex = profileIndex
1349 1349
1350 1350 dataOut.ippSeconds /= self.__nTxs
1351 1351
1352 1352 return dataOut
1353 1353
1354 1354 class SplitProfiles(Operation):
1355 1355
1356 1356 def __init__(self, **kwargs):
1357 1357
1358 1358 Operation.__init__(self, **kwargs)
1359 1359
1360 1360 def run(self, dataOut, n):
1361 1361
1362 1362 dataOut.flagNoData = True
1363 1363 profileIndex = None
1364 1364
1365 1365 if dataOut.flagDataAsBlock:
1366 1366
1367 1367 #nchannels, nprofiles, nsamples
1368 1368 shape = dataOut.data.shape
1369 1369
1370 1370 if shape[2] % n != 0:
1371 1371 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1372 1372
1373 1373 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1374 1374
1375 1375 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1376 1376 dataOut.flagNoData = False
1377 1377
1378 1378 profileIndex = int(dataOut.nProfiles/n) - 1
1379 1379
1380 1380 else:
1381 1381
1382 1382 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1383 1383
1384 1384 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1385 1385
1386 1386 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1387 1387
1388 1388 dataOut.nProfiles = int(dataOut.nProfiles*n)
1389 1389
1390 1390 dataOut.profileIndex = profileIndex
1391 1391
1392 1392 dataOut.ippSeconds /= n
1393 1393
1394 1394 return dataOut
1395 1395
1396 1396 class CombineProfiles(Operation):
1397 1397 def __init__(self, **kwargs):
1398 1398
1399 1399 Operation.__init__(self, **kwargs)
1400 1400
1401 1401 self.__remData = None
1402 1402 self.__profileIndex = 0
1403 1403
1404 1404 def run(self, dataOut, n):
1405 1405
1406 1406 dataOut.flagNoData = True
1407 1407 profileIndex = None
1408 1408
1409 1409 if dataOut.flagDataAsBlock:
1410 1410
1411 1411 #nchannels, nprofiles, nsamples
1412 1412 shape = dataOut.data.shape
1413 1413 new_shape = shape[0], shape[1]/n, shape[2]*n
1414 1414
1415 1415 if shape[1] % n != 0:
1416 1416 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1417 1417
1418 1418 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1419 1419 dataOut.flagNoData = False
1420 1420
1421 1421 profileIndex = int(dataOut.nProfiles*n) - 1
1422 1422
1423 1423 else:
1424 1424
1425 1425 #nchannels, nsamples
1426 1426 if self.__remData is None:
1427 1427 newData = dataOut.data
1428 1428 else:
1429 1429 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1430 1430
1431 1431 self.__profileIndex += 1
1432 1432
1433 1433 if self.__profileIndex < n:
1434 1434 self.__remData = newData
1435 1435 #continue
1436 1436 return
1437 1437
1438 1438 self.__profileIndex = 0
1439 1439 self.__remData = None
1440 1440
1441 1441 dataOut.data = newData
1442 1442 dataOut.flagNoData = False
1443 1443
1444 1444 profileIndex = dataOut.profileIndex/n
1445 1445
1446 1446
1447 1447 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1448 1448
1449 1449 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1450 1450
1451 1451 dataOut.nProfiles = int(dataOut.nProfiles/n)
1452 1452
1453 1453 dataOut.profileIndex = profileIndex
1454 1454
1455 1455 dataOut.ippSeconds *= n
1456 1456
1457 1457 return dataOut
1458 1458
1459 1459 class PulsePairVoltage(Operation):
1460 1460 '''
1461 1461 Function PulsePair(Signal Power, Velocity)
1462 1462 The real component of Lag[0] provides Intensity Information
1463 1463 The imag component of Lag[1] Phase provides Velocity Information
1464 1464
1465 1465 Configuration Parameters:
1466 1466 nPRF = Number of Several PRF
1467 1467 theta = Degree Azimuth angel Boundaries
1468 1468
1469 1469 Input:
1470 1470 self.dataOut
1471 1471 lag[N]
1472 1472 Affected:
1473 1473 self.dataOut.spc
1474 1474 '''
1475 1475 isConfig = False
1476 1476 __profIndex = 0
1477 1477 __initime = None
1478 1478 __lastdatatime = None
1479 1479 __buffer = None
1480 1480 noise = None
1481 1481 __dataReady = False
1482 1482 n = None
1483 1483 __nch = 0
1484 1484 __nHeis = 0
1485 1485 removeDC = False
1486 1486 ipp = None
1487 1487 lambda_ = 0
1488 1488
1489 1489 def __init__(self,**kwargs):
1490 1490 Operation.__init__(self,**kwargs)
1491 1491
1492 1492 def setup(self, dataOut, n = None, removeDC=False):
1493 1493 '''
1494 1494 n= Numero de PRF's de entrada
1495 1495 '''
1496 1496 self.__initime = None
1497 1497 self.__lastdatatime = 0
1498 1498 self.__dataReady = False
1499 1499 self.__buffer = 0
1500 1500 self.__profIndex = 0
1501 1501 self.noise = None
1502 1502 self.__nch = dataOut.nChannels
1503 1503 self.__nHeis = dataOut.nHeights
1504 1504 self.removeDC = removeDC
1505 1505 self.lambda_ = 3.0e8/(9345.0e6)
1506 1506 self.ippSec = dataOut.ippSeconds
1507 1507 self.nCohInt = dataOut.nCohInt
1508 1508
1509 1509 if n == None:
1510 1510 raise ValueError("n should be specified.")
1511 1511
1512 1512 if n != None:
1513 1513 if n<2:
1514 1514 raise ValueError("n should be greater than 2")
1515 1515
1516 1516 self.n = n
1517 1517 self.__nProf = n
1518 1518
1519 1519 self.__buffer = numpy.zeros((dataOut.nChannels,
1520 1520 n,
1521 1521 dataOut.nHeights),
1522 1522 dtype='complex')
1523 1523
1524 1524 def putData(self,data):
1525 1525 '''
1526 1526 Add a profile to he __buffer and increase in one the __profiel Index
1527 1527 '''
1528 1528 self.__buffer[:,self.__profIndex,:]= data
1529 1529 self.__profIndex += 1
1530 1530 return
1531 1531
1532 1532 def pushData(self,dataOut):
1533 1533 '''
1534 1534 Return the PULSEPAIR and the profiles used in the operation
1535 1535 Affected : self.__profileIndex
1536 1536 '''
1537 1537 #----------------- Remove DC-----------------------------------
1538 1538 if self.removeDC==True:
1539 1539 mean = numpy.mean(self.__buffer,1)
1540 1540 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1541 1541 dc= numpy.tile(tmp,[1,self.__nProf,1])
1542 1542 self.__buffer = self.__buffer - dc
1543 1543 #------------------Calculo de Potencia ------------------------
1544 1544 pair0 = self.__buffer*numpy.conj(self.__buffer)
1545 1545 pair0 = pair0.real
1546 1546 lag_0 = numpy.sum(pair0,1)
1547 1547 #------------------Calculo de Ruido x canal--------------------
1548 1548 self.noise = numpy.zeros(self.__nch)
1549 1549 for i in range(self.__nch):
1550 1550 daux = numpy.sort(pair0[i,:,:],axis= None)
1551 1551 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1552 1552
1553 1553 self.noise = self.noise.reshape(self.__nch,1)
1554 1554 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1555 1555 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1556 1556 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1557 1557 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1558 1558 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1559 1559 #-------------------- Power --------------------------------------------------
1560 1560 data_power = lag_0/(self.n*self.nCohInt)
1561 1561 #------------------ Senal ---------------------------------------------------
1562 1562 data_intensity = pair0 - noise_buffer
1563 1563 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1564 1564 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1565 1565 for i in range(self.__nch):
1566 1566 for j in range(self.__nHeis):
1567 1567 if data_intensity[i][j] < 0:
1568 1568 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1569 1569
1570 1570 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1571 1571 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1572 1572 lag_1 = numpy.sum(pair1,1)
1573 1573 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1574 1574 data_velocity = (self.lambda_/2.0)*data_freq
1575 1575
1576 1576 #---------------- Potencia promedio estimada de la Senal-----------
1577 1577 lag_0 = lag_0/self.n
1578 1578 S = lag_0-self.noise
1579 1579
1580 1580 #---------------- Frecuencia Doppler promedio ---------------------
1581 1581 lag_1 = lag_1/(self.n-1)
1582 1582 R1 = numpy.abs(lag_1)
1583 1583
1584 1584 #---------------- Calculo del SNR----------------------------------
1585 1585 data_snrPP = S/self.noise
1586 1586 for i in range(self.__nch):
1587 1587 for j in range(self.__nHeis):
1588 1588 if data_snrPP[i][j] < 1.e-20:
1589 1589 data_snrPP[i][j] = 1.e-20
1590 1590
1591 1591 #----------------- Calculo del ancho espectral ----------------------
1592 1592 L = S/R1
1593 1593 L = numpy.where(L<0,1,L)
1594 1594 L = numpy.log(L)
1595 1595 tmp = numpy.sqrt(numpy.absolute(L))
1596 1596 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1597 1597 n = self.__profIndex
1598 1598
1599 1599 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1600 1600 self.__profIndex = 0
1601 1601 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1602 1602
1603 1603
1604 1604 def pulsePairbyProfiles(self,dataOut):
1605 1605
1606 1606 self.__dataReady = False
1607 1607 data_power = None
1608 1608 data_intensity = None
1609 1609 data_velocity = None
1610 1610 data_specwidth = None
1611 1611 data_snrPP = None
1612 1612 self.putData(data=dataOut.data)
1613 1613 if self.__profIndex == self.n:
1614 1614 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1615 1615 self.__dataReady = True
1616 1616
1617 1617 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1618 1618
1619 1619
1620 1620 def pulsePairOp(self, dataOut, datatime= None):
1621 1621
1622 1622 if self.__initime == None:
1623 1623 self.__initime = datatime
1624 1624 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1625 1625 self.__lastdatatime = datatime
1626 1626
1627 1627 if data_power is None:
1628 1628 return None, None, None,None,None,None
1629 1629
1630 1630 avgdatatime = self.__initime
1631 1631 deltatime = datatime - self.__lastdatatime
1632 1632 self.__initime = datatime
1633 1633
1634 1634 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1635 1635
1636 1636 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1637 1637
1638 1638 if not self.isConfig:
1639 1639 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1640 1640 self.isConfig = True
1641 1641 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1642 1642 dataOut.flagNoData = True
1643 1643
1644 1644 if self.__dataReady:
1645 1645 dataOut.nCohInt *= self.n
1646 1646 dataOut.dataPP_POW = data_intensity # S
1647 1647 dataOut.dataPP_POWER = data_power # P
1648 1648 dataOut.dataPP_DOP = data_velocity
1649 1649 dataOut.dataPP_SNR = data_snrPP
1650 1650 dataOut.dataPP_WIDTH = data_specwidth
1651 1651 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1652 1652 dataOut.utctime = avgdatatime
1653 1653 dataOut.flagNoData = False
1654 1654 return dataOut
1655 1655
1656 1656
1657 1657
1658 1658 # import collections
1659 1659 # from scipy.stats import mode
1660 1660 #
1661 1661 # class Synchronize(Operation):
1662 1662 #
1663 1663 # isConfig = False
1664 1664 # __profIndex = 0
1665 1665 #
1666 1666 # def __init__(self, **kwargs):
1667 1667 #
1668 1668 # Operation.__init__(self, **kwargs)
1669 1669 # # self.isConfig = False
1670 1670 # self.__powBuffer = None
1671 1671 # self.__startIndex = 0
1672 1672 # self.__pulseFound = False
1673 1673 #
1674 1674 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1675 1675 #
1676 1676 # #Read data
1677 1677 #
1678 1678 # powerdB = dataOut.getPower(channel = channel)
1679 1679 # noisedB = dataOut.getNoise(channel = channel)[0]
1680 1680 #
1681 1681 # self.__powBuffer.extend(powerdB.flatten())
1682 1682 #
1683 1683 # dataArray = numpy.array(self.__powBuffer)
1684 1684 #
1685 1685 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1686 1686 #
1687 1687 # maxValue = numpy.nanmax(filteredPower)
1688 1688 #
1689 1689 # if maxValue < noisedB + 10:
1690 1690 # #No se encuentra ningun pulso de transmision
1691 1691 # return None
1692 1692 #
1693 1693 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1694 1694 #
1695 1695 # if len(maxValuesIndex) < 2:
1696 1696 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1697 1697 # return None
1698 1698 #
1699 1699 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1700 1700 #
1701 1701 # #Seleccionar solo valores con un espaciamiento de nSamples
1702 1702 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1703 1703 #
1704 1704 # if len(pulseIndex) < 2:
1705 1705 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1706 1706 # return None
1707 1707 #
1708 1708 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1709 1709 #
1710 1710 # #remover senales que se distancien menos de 10 unidades o muestras
1711 1711 # #(No deberian existir IPP menor a 10 unidades)
1712 1712 #
1713 1713 # realIndex = numpy.where(spacing > 10 )[0]
1714 1714 #
1715 1715 # if len(realIndex) < 2:
1716 1716 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1717 1717 # return None
1718 1718 #
1719 1719 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1720 1720 # realPulseIndex = pulseIndex[realIndex]
1721 1721 #
1722 1722 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1723 1723 #
1724 1724 # print "IPP = %d samples" %period
1725 1725 #
1726 1726 # self.__newNSamples = dataOut.nHeights #int(period)
1727 1727 # self.__startIndex = int(realPulseIndex[0])
1728 1728 #
1729 1729 # return 1
1730 1730 #
1731 1731 #
1732 1732 # def setup(self, nSamples, nChannels, buffer_size = 4):
1733 1733 #
1734 1734 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1735 1735 # maxlen = buffer_size*nSamples)
1736 1736 #
1737 1737 # bufferList = []
1738 1738 #
1739 1739 # for i in range(nChannels):
1740 1740 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1741 1741 # maxlen = buffer_size*nSamples)
1742 1742 #
1743 1743 # bufferList.append(bufferByChannel)
1744 1744 #
1745 1745 # self.__nSamples = nSamples
1746 1746 # self.__nChannels = nChannels
1747 1747 # self.__bufferList = bufferList
1748 1748 #
1749 1749 # def run(self, dataOut, channel = 0):
1750 1750 #
1751 1751 # if not self.isConfig:
1752 1752 # nSamples = dataOut.nHeights
1753 1753 # nChannels = dataOut.nChannels
1754 1754 # self.setup(nSamples, nChannels)
1755 1755 # self.isConfig = True
1756 1756 #
1757 1757 # #Append new data to internal buffer
1758 1758 # for thisChannel in range(self.__nChannels):
1759 1759 # bufferByChannel = self.__bufferList[thisChannel]
1760 1760 # bufferByChannel.extend(dataOut.data[thisChannel])
1761 1761 #
1762 1762 # if self.__pulseFound:
1763 1763 # self.__startIndex -= self.__nSamples
1764 1764 #
1765 1765 # #Finding Tx Pulse
1766 1766 # if not self.__pulseFound:
1767 1767 # indexFound = self.__findTxPulse(dataOut, channel)
1768 1768 #
1769 1769 # if indexFound == None:
1770 1770 # dataOut.flagNoData = True
1771 1771 # return
1772 1772 #
1773 1773 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1774 1774 # self.__pulseFound = True
1775 1775 # self.__startIndex = indexFound
1776 1776 #
1777 1777 # #If pulse was found ...
1778 1778 # for thisChannel in range(self.__nChannels):
1779 1779 # bufferByChannel = self.__bufferList[thisChannel]
1780 1780 # #print self.__startIndex
1781 1781 # x = numpy.array(bufferByChannel)
1782 1782 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1783 1783 #
1784 1784 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1785 1785 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1786 1786 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1787 1787 #
1788 1788 # dataOut.data = self.__arrayBuffer
1789 1789 #
1790 1790 # self.__startIndex += self.__newNSamples
1791 1791 #
1792 1792 # return
1793 1793 class SSheightProfiles(Operation):
1794 1794
1795 1795 step = None
1796 1796 nsamples = None
1797 1797 bufferShape = None
1798 1798 profileShape = None
1799 1799 sshProfiles = None
1800 1800 profileIndex = None
1801 1801
1802 1802 def __init__(self, **kwargs):
1803 1803
1804 1804 Operation.__init__(self, **kwargs)
1805 1805 self.isConfig = False
1806 1806
1807 1807 def setup(self,dataOut ,step = None , nsamples = None):
1808 1808
1809 1809 if step == None and nsamples == None:
1810 1810 raise ValueError("step or nheights should be specified ...")
1811 1811
1812 1812 self.step = step
1813 1813 self.nsamples = nsamples
1814 1814 self.__nChannels = dataOut.nChannels
1815 1815 self.__nProfiles = dataOut.nProfiles
1816 1816 self.__nHeis = dataOut.nHeights
1817 1817 shape = dataOut.data.shape #nchannels, nprofiles, nsamples
1818 1818
1819 1819 residue = (shape[1] - self.nsamples) % self.step
1820 1820 if residue != 0:
1821 1821 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
1822 1822
1823 1823 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1824 1824 numberProfile = self.nsamples
1825 1825 numberSamples = (shape[1] - self.nsamples)/self.step
1826 1826
1827 1827 self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles
1828 1828 self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples
1829 1829
1830 1830 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1831 1831 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1832 1832
1833 1833 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1834 1834 dataOut.flagNoData = True
1835 1835
1836 1836 profileIndex = None
1837 1837 #print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights)
1838 1838 #print(dataOut.getFreqRange(1)/1000.)
1839 1839 #exit(1)
1840 1840 if dataOut.flagDataAsBlock:
1841 1841 dataOut.data = numpy.average(dataOut.data,axis=1)
1842 1842 #print("jee")
1843 1843 dataOut.flagDataAsBlock = False
1844 1844 if not self.isConfig:
1845 1845 self.setup(dataOut, step=step , nsamples=nsamples)
1846 1846 #print("Setup done")
1847 1847 self.isConfig = True
1848 1848
1849 1849
1850 1850 if code is not None:
1851 1851 code = numpy.array(code)
1852 1852 code_block = code
1853 1853
1854 1854 if repeat is not None:
1855 1855 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1856 1856 #print(code_block.shape)
1857 1857 for i in range(self.buffer.shape[1]):
1858 1858
1859 1859 if code is not None:
1860 1860 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block
1861 1861
1862 1862 else:
1863 1863
1864 1864 self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1865 1865
1866 1866 #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
1867 1867
1868 1868 for j in range(self.buffer.shape[0]):
1869 1869 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1870 1870
1871 1871 profileIndex = self.nsamples
1872 1872 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1873 1873 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1874 1874 #print("ippSeconds, dH: ",ippSeconds,deltaHeight)
1875 1875 try:
1876 1876 if dataOut.concat_m is not None:
1877 1877 ippSeconds= ippSeconds/float(dataOut.concat_m)
1878 1878 #print "Profile concat %d"%dataOut.concat_m
1879 1879 except:
1880 1880 pass
1881 1881
1882 1882 dataOut.data = self.sshProfiles
1883 1883 dataOut.flagNoData = False
1884 1884 dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
1885 1885 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1886 1886
1887 1887 dataOut.profileIndex = profileIndex
1888 1888 dataOut.flagDataAsBlock = True
1889 1889 dataOut.ippSeconds = ippSeconds
1890 1890 dataOut.step = self.step
1891 1891 #print(numpy.shape(dataOut.data))
1892 1892 #exit(1)
1893 1893 #print("new data shape and time:", dataOut.data.shape, dataOut.utctime)
1894 1894
1895 1895 return dataOut
1896 1896 ################################################################################3############################3
1897 1897 ################################################################################3############################3
1898 1898 ################################################################################3############################3
1899 1899 ################################################################################3############################3
1900 1900
1901 1901 class SSheightProfiles2(Operation):
1902 1902 '''
1903 1903 Procesa por perfiles y por bloques
1904 1904 VersiΓ³n corregida y actualizada para trabajar con RemoveProfileSats2
1905 1905 Usar esto
1906 1906 '''
1907 1907
1908 1908
1909 1909 bufferShape = None
1910 1910 profileShape = None
1911 1911 sshProfiles = None
1912 1912 profileIndex = None
1913 1913 #nsamples = None
1914 1914 #step = None
1915 1915 #deltaHeight = None
1916 1916 #init_range = None
1917 1917 __slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels',
1918 1918 '__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights')
1919 1919
1920 1920 def __init__(self, **kwargs):
1921 1921
1922 1922 Operation.__init__(self, **kwargs)
1923 1923 self.isConfig = False
1924 1924
1925 1925 def setup(self,dataOut ,step = None , nsamples = None):
1926 1926
1927 1927 if step == None and nsamples == None:
1928 1928 raise ValueError("step or nheights should be specified ...")
1929 1929
1930 1930 self.step = step
1931 1931 self.nsamples = nsamples
1932 1932 self.__nChannels = int(dataOut.nChannels)
1933 1933 self.__nProfiles = int(dataOut.nProfiles)
1934 1934 self.__nHeis = int(dataOut.nHeights)
1935 1935
1936 1936 residue = (self.__nHeis - self.nsamples) % self.step
1937 1937 if residue != 0:
1938 1938 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,self.__nProfiles - self.nsamples,residue))
1939 1939
1940 1940 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1941 1941 self.init_range = dataOut.heightList[0]
1942 1942 #numberProfile = self.nsamples
1943 1943 numberSamples = (self.__nHeis - self.nsamples)/self.step
1944 1944
1945 1945 self.new_nHeights = numberSamples
1946 1946
1947 1947 self.bufferShape = int(self.__nChannels), int(numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
1948 1948 self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples) # nchannels, nprofiles, nsamples
1949 1949
1950 1950 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1951 1951 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1952 1952
1953 1953 def getNewProfiles(self, data, code=None, repeat=None):
1954 1954
1955 1955 if code is not None:
1956 1956 code = numpy.array(code)
1957 1957 code_block = code
1958 1958
1959 1959 if repeat is not None:
1960 1960 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1961 1961 if data.ndim < 3:
1962 1962 data = data.reshape(self.__nChannels,1,self.__nHeis )
1963 1963 #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape, code_block.shape)
1964 1964 for ch in range(self.__nChannels):
1965 1965 for i in range(int(self.new_nHeights)): #nuevas alturas
1966 1966 if code is not None:
1967 1967 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]*code_block
1968 1968 else:
1969 1969 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1970 1970
1971 1971 for j in range(self.__nChannels): #en los cananles
1972 1972 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
1973 1973 #print("new profs Done")
1974 1974
1975 1975
1976 1976
1977 1977 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1978 1978 # print("running")
1979 1979 if dataOut.flagNoData == True:
1980 1980 return dataOut
1981 1981 dataOut.flagNoData = True
1982 1982 #print("init data shape:", dataOut.data.shape)
1983 1983 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
1984 1984 # int(dataOut.nProfiles),int(dataOut.nHeights)))
1985 1985
1986 1986 profileIndex = None
1987 1987 # if not dataOut.flagDataAsBlock:
1988 1988 # dataOut.nProfiles = 1
1989 1989
1990 1990 if not self.isConfig:
1991 1991 self.setup(dataOut, step=step , nsamples=nsamples)
1992 1992 #print("Setup done")
1993 1993 self.isConfig = True
1994 1994
1995 1995 dataBlock = None
1996 1996
1997 1997 nprof = 1
1998 1998 if dataOut.flagDataAsBlock:
1999 1999 nprof = int(dataOut.nProfiles)
2000 2000
2001 2001 #print("dataOut nProfiles:", dataOut.nProfiles)
2002 2002 for profile in range(nprof):
2003 2003 if dataOut.flagDataAsBlock:
2004 2004 #print("read blocks")
2005 2005 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
2006 2006 else:
2007 2007 #print("read profiles")
2008 2008 self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe
2009 2009 if profile == 0:
2010 2010 dataBlock = self.sshProfiles.copy()
2011 2011 else: #by blocks
2012 2012 dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
2013 2013 #print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
2014 2014
2015 2015 profileIndex = self.nsamples
2016 2016 #deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
2017 2017 ippSeconds = (self.deltaHeight*1.0e-6)/(0.15)
2018 2018
2019 2019
2020 2020 dataOut.data = dataBlock
2021 2021 #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
2022 2022 dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
2023 2023 dataOut.sampled_heightsFFT = self.nsamples
2024 2024 dataOut.ippSeconds = ippSeconds
2025 2025 dataOut.step = self.step
2026 2026 dataOut.deltaHeight = self.step*self.deltaHeight
2027 2027 dataOut.flagNoData = False
2028 2028 if dataOut.flagDataAsBlock:
2029 2029 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
2030 2030
2031 2031 else:
2032 2032 dataOut.nProfiles = int(self.nsamples)
2033 2033 dataOut.profileIndex = dataOut.nProfiles
2034 2034 dataOut.flagDataAsBlock = True
2035 2035
2036 2036 dataBlock = None
2037 2037
2038 2038 #print("new data shape:", dataOut.data.shape, dataOut.utctime)
2039 2039
2040 2040 #update Processing Header:
2041 2041 dataOut.processingHeaderObj.heightList = dataOut.heightList
2042 2042 dataOut.processingHeaderObj.ipp = ippSeconds
2043 2043 dataOut.processingHeaderObj.heightResolution = dataOut.deltaHeight
2044 2044 #dataOut.processingHeaderObj.profilesPerBlock = nProfiles
2045 2045
2046 2046 # # dataOut.data = CH, PROFILES, HEIGHTS
2047 2047 #print(dataOut.data .shape)
2048 2048 if dataOut.flagProfilesByRange:
2049 2049 # #assuming the same remotion for all channels
2050 2050 aux = [ self.nsamples - numpy.count_nonzero(dataOut.data[0, :, h]==0) for h in range(len(dataOut.heightList))]
2051 2051 dataOut.nProfilesByRange = (numpy.asarray(aux)).reshape((1,len(dataOut.heightList) ))
2052 2052 #print(dataOut.nProfilesByRange.shape)
2053 2053 else:
2054 2054 dataOut.nProfilesByRange = numpy.ones((1, len(dataOut.heightList)))*dataOut.nProfiles
2055 2055 return dataOut
2056 2056
2057 2057
2058 2058
2059 2059 class RemoveProfileSats(Operation):
2060 2060 '''
2061 2061 Escrito: Joab Apaza
2062 2062
2063 2063 Omite los perfiles contaminados con seΓ±al de satΓ©lites, usando una altura de referencia
2064 2064 In: minHei = min_sat_range
2065 2065 max_sat_range
2066 2066 min_hei_ref
2067 2067 max_hei_ref
2068 2068 th = diference between profiles mean, ref and sats
2069 2069 Out:
2070 2070 profile clean
2071 2071 '''
2072 2072
2073 2073
2074 2074 __buffer_data = []
2075 2075 __buffer_times = []
2076 2076
2077 2077 buffer = None
2078 2078
2079 2079 outliers_IDs_list = []
2080 2080
2081 2081
2082 2082 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
2083 2083 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
2084 2084 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'thdB')
2085 2085 def __init__(self, **kwargs):
2086 2086
2087 2087 Operation.__init__(self, **kwargs)
2088 2088 self.isConfig = False
2089 2089
2090 2090 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=15,
2091 2091 minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
2092 2092
2093 2093 if n == None and timeInterval == None:
2094 2094 raise ValueError("nprofiles or timeInterval should be specified ...")
2095 2095
2096 2096 if n != None:
2097 2097 self.n = n
2098 2098
2099 2099 self.navg = navg
2100 2100 self.profileMargin = profileMargin
2101 2101 self.thHistOutlier = thHistOutlier
2102 2102 self.__profIndex = 0
2103 2103 self.buffer = None
2104 2104 self._ipp = dataOut.ippSeconds
2105 2105 self.n_prof_released = 0
2106 2106 self.heightList = dataOut.heightList
2107 2107 self.init_prof = 0
2108 2108 self.end_prof = 0
2109 2109 self.__count_exec = 0
2110 2110 self.__profIndex = 0
2111 2111 self.first_utcBlock = None
2112 2112 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
2113 2113 minHei = minHei
2114 2114 maxHei = maxHei
2115 2115 if minHei==None :
2116 2116 minHei = dataOut.heightList[0]
2117 2117 if maxHei==None :
2118 2118 maxHei = dataOut.heightList[-1]
2119 2119 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
2120 2120 self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
2121 2121 self.nChannels = dataOut.nChannels
2122 2122 self.nHeights = dataOut.nHeights
2123 2123 self.test_counter = 0
2124 2124 self.thdB = thdB
2125 2125
2126 2126 def filterSatsProfiles(self):
2127 2127 data = self.__buffer_data
2128 2128 #print(data.shape)
2129 2129 nChannels, profiles, heights = data.shape
2130 2130 indexes=numpy.zeros([], dtype=int)
2131 2131 outliers_IDs=[]
2132 2132 for c in range(nChannels):
2133 2133 #print(self.min_ref,self.max_ref)
2134 2134 noise_ref = 10* numpy.log10((data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])).real)
2135 2135 #print("Noise ",numpy.percentile(noise_ref,95))
2136 2136 p95 = numpy.percentile(noise_ref,95)
2137 2137 noise_ref = noise_ref.mean()
2138 2138 #print("Noise ",noise_ref
2139 2139
2140 2140
2141 2141 for h in range(self.minHei_idx, self.maxHei_idx):
2142 2142 power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
2143 2143 #th = noise_ref + self.thdB
2144 2144 th = noise_ref + 1.5*(p95-noise_ref)
2145 2145 index = numpy.where(power > th )
2146 2146 if index[0].size > 10 and index[0].size < int(self.navg*profiles):
2147 2147 indexes = numpy.append(indexes, index[0])
2148 2148 #print(index[0])
2149 2149 #print(index[0])
2150 2150
2151 2151 # fig,ax = plt.subplots()
2152 2152 # #ax.set_title(str(k)+" "+str(j))
2153 2153 # x=range(len(power))
2154 2154 # ax.scatter(x,power)
2155 2155 # #ax.axvline(index)
2156 2156 # plt.grid()
2157 2157 # plt.show()
2158 2158 #print(indexes)
2159 2159
2160 2160 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2161 2161 #outliers_IDs = numpy.unique(outliers_IDs)
2162 2162
2163 2163 outs_lines = numpy.unique(indexes)
2164 2164
2165 2165
2166 2166 #Agrupando el histograma de outliers,
2167 2167 my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=True)
2168 2168
2169 2169
2170 2170 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
2171 2171 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
2172 2172 hist_outliers_indexes = hist_outliers_indexes[0]
2173 2173 # if len(hist_outliers_indexes>0):
2174 2174 # hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
2175 2175 #print(hist_outliers_indexes)
2176 2176 #print(bins, hist_outliers_indexes)
2177 2177 bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
2178 2178 outlier_loc_index = []
2179 2179 # for n in range(len(bins_outliers_indexes)):
2180 2180 # for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n]+ self.profileMargin):
2181 2181 # outlier_loc_index.append(e)
2182 2182 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n]+ profiles//100 + self.profileMargin) ]
2183 2183 outlier_loc_index = numpy.asarray(outlier_loc_index)
2184 2184
2185 2185
2186 2186
2187 2187
2188 2188 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
2189 2189 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
2190 2190 #print("outliers final: ", outlier_loc_index)
2191 2191
2192 2192 from matplotlib import pyplot as plt
2193 2193 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
2194 2194 fig, ax = plt.subplots(1,2,figsize=(8, 6))
2195 2195 dat = data[0,:,:].real
2196 2196 dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
2197 2197 m = numpy.nanmean(dat)
2198 2198 o = numpy.nanstd(dat)
2199 2199 #print(m, o, x.shape, y.shape)
2200 2200 #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2201 2201 c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = 50, vmax = 75)
2202 2202 ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
2203 2203 fig.colorbar(c)
2204 2204 ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
2205 2205 ax[1].hist(outs_lines,bins=my_bins)
2206 2206 plt.show()
2207 2207
2208 2208
2209 2209 self.outliers_IDs_list = outlier_loc_index
2210 2210 #print("outs list: ", self.outliers_IDs_list)
2211 2211 return data
2212 2212
2213 2213
2214 2214
2215 2215 def fillBuffer(self, data, datatime):
2216 2216
2217 2217 if self.__profIndex == 0:
2218 2218 self.__buffer_data = data.copy()
2219 2219
2220 2220 else:
2221 2221 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2222 2222 self.__profIndex += 1
2223 2223 self.__buffer_times.append(datatime)
2224 2224
2225 2225 def getData(self, data, datatime=None):
2226 2226
2227 2227 if self.__profIndex == 0:
2228 2228 self.__initime = datatime
2229 2229
2230 2230
2231 2231 self.__dataReady = False
2232 2232
2233 2233 self.fillBuffer(data, datatime)
2234 2234 dataBlock = None
2235 2235
2236 2236 if self.__profIndex == self.n:
2237 2237 #print("apnd : ",data)
2238 2238 dataBlock = self.filterSatsProfiles()
2239 2239 self.__dataReady = True
2240 2240
2241 2241 return dataBlock
2242 2242
2243 2243 if dataBlock is None:
2244 2244 return None, None
2245 2245
2246 2246
2247 2247
2248 2248 return dataBlock
2249 2249
2250 2250 def releaseBlock(self):
2251 2251
2252 2252 if self.n % self.lenProfileOut != 0:
2253 2253 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2254 2254 return None
2255 2255
2256 2256 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2257 2257
2258 2258 self.init_prof = self.end_prof
2259 2259 self.end_prof += self.lenProfileOut
2260 2260 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2261 2261 self.n_prof_released += 1
2262 2262
2263 2263 return data
2264 2264
2265 2265 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,
2266 2266 th_hist_outlier=15,minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
2267 2267
2268 2268 if not self.isConfig:
2269 2269 #print("init p idx: ", dataOut.profileIndex )
2270 2270 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,
2271 2271 minHei=minHei, maxHei=maxHei, minRef=minRef, maxRef=maxRef, thdB=thdB)
2272 2272 self.isConfig = True
2273 2273
2274 2274 dataBlock = None
2275 2275
2276 2276 if not dataOut.buffer_empty: #hay datos acumulados
2277 2277
2278 2278 if self.init_prof == 0:
2279 2279 self.n_prof_released = 0
2280 2280 self.lenProfileOut = nProfilesOut
2281 2281 dataOut.flagNoData = False
2282 2282 #print("tp 2 ",dataOut.data.shape)
2283 2283
2284 2284 self.init_prof = 0
2285 2285 self.end_prof = self.lenProfileOut
2286 2286
2287 2287 dataOut.nProfiles = self.lenProfileOut
2288 2288 if nProfilesOut == 1:
2289 2289 dataOut.flagDataAsBlock = False
2290 2290 else:
2291 2291 dataOut.flagDataAsBlock = True
2292 2292 #print("prof: ",self.init_prof)
2293 2293 dataOut.flagNoData = False
2294 2294 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
2295 2295 #print("omitting: ", self.n_prof_released)
2296 2296 dataOut.flagNoData = True
2297 2297 dataOut.ippSeconds = self._ipp
2298 2298 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
2299 2299 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
2300 2300 #dataOut.data = self.releaseBlock()
2301 2301 #########################################################3
2302 2302 if self.n % self.lenProfileOut != 0:
2303 2303 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2304 2304 return None
2305 2305
2306 2306 dataOut.data = None
2307 2307
2308 2308 if nProfilesOut == 1:
2309 2309 dataOut.data = self.buffer[:,self.end_prof-1,:] #ch, prof, alt
2310 2310 else:
2311 2311 dataOut.data = self.buffer[:,self.init_prof:self.end_prof,:] #ch, prof, alt
2312 2312
2313 2313 self.init_prof = self.end_prof
2314 2314 self.end_prof += self.lenProfileOut
2315 2315 #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
2316 2316 self.n_prof_released += 1
2317 2317
2318 2318 if self.end_prof >= (self.n +self.lenProfileOut):
2319 2319
2320 2320 self.init_prof = 0
2321 2321 self.__profIndex = 0
2322 2322 self.buffer = None
2323 2323 dataOut.buffer_empty = True
2324 2324 self.outliers_IDs_list = []
2325 2325 self.n_prof_released = 0
2326 2326 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
2327 2327 #print("cleaning...", dataOut.buffer_empty)
2328 2328 dataOut.profileIndex = 0 #self.lenProfileOut
2329 2329 ####################################################################
2330 2330 return dataOut
2331 2331
2332 2332
2333 2333 #print("tp 223 ",dataOut.data.shape)
2334 2334 dataOut.flagNoData = True
2335 2335
2336 2336
2337 2337
2338 2338 try:
2339 2339 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
2340 2340 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
2341 2341 self.__count_exec +=1
2342 2342 except Exception as e:
2343 2343 print("Error getting profiles data",self.__count_exec )
2344 2344 print(e)
2345 2345 sys.exit()
2346 2346
2347 2347 if self.__dataReady:
2348 2348 #print("omitting: ", len(self.outliers_IDs_list))
2349 2349 self.__count_exec = 0
2350 2350 #dataOut.data =
2351 2351 #self.buffer = numpy.flip(dataBlock, axis=1)
2352 2352 self.buffer = dataBlock
2353 2353 self.first_utcBlock = self.__initime
2354 2354 dataOut.utctime = self.__initime
2355 2355 dataOut.nProfiles = self.__profIndex
2356 2356 #dataOut.flagNoData = False
2357 2357 self.init_prof = 0
2358 2358 self.__profIndex = 0
2359 2359 self.__initime = None
2360 2360 dataBlock = None
2361 2361 self.__buffer_times = []
2362 2362 dataOut.error = False
2363 2363 dataOut.useInputBuffer = True
2364 2364 dataOut.buffer_empty = False
2365 2365 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
2366 2366
2367 2367
2368 2368
2369 2369 #print(self.__count_exec)
2370 2370
2371 2371 return dataOut
2372 2372
2373 2373
2374 2374 class RemoveProfileSats2(Operation):
2375 2375 '''
2376 2376 Escrito: Joab Apaza
2377 2377
2378 2378 Omite los perfiles contaminados con seΓ±al de satΓ©lites, usando una altura de referencia
2379 2379 promedia todas las alturas para los cΓ‘lculos
2380 2380 In:
2381 2381 n = Cantidad de perfiles que se acumularan, usualmente 10 segundos
2382 2382 navg = Porcentaje de perfiles que puede considerarse como satΓ©lite, mΓ‘ximo 90%
2383 2383 minHei =
2384 2384 minRef =
2385 2385 maxRef =
2386 2386 nBins =
2387 2387 profile_margin =
2388 2388 th_hist_outlier =
2389 2389 nProfilesOut =
2390 2390
2391 2391 Pensado para remover interferencias de las YAGI, se puede adaptar a otras interferencias
2392 2392
2393 2393 remYagi = Activa la funcion de remociΓ³n de interferencias de la YAGI
2394 2394 nProfYagi = Cantidad de perfiles que son afectados, acorde NTX de la YAGI
2395 2395 offYagi =
2396 2396 minHJULIA = Altura mΓ­nima donde aparece la seΓ±al referencia de JULIA (-50)
2397 2397 maxHJULIA = Altura mΓ‘xima donde aparece la seΓ±al referencia de JULIA (-15)
2398 2398
2399 2399 debug = Activa los grΓ‘ficos, recomendable ejecutar para ajustar los parΓ‘metros
2400 2400 para un experimento en especΓ­fico.
2401 2401
2402 2402 ** se modifica para remover interferencias puntuales, es decir, desde otros radares.
2403 2403 Inicialmente se ha configurado para omitir tambiΓ©n los perfiles de la YAGI en los datos
2404 2404 de AMISR-ISR.
2405 2405
2406 2406 Out:
2407 2407 profile clean
2408 2408 '''
2409 2409
2410 2410
2411 2411 __buffer_data = []
2412 2412 __buffer_times = []
2413 2413
2414 2414 buffer = None
2415 2415
2416 2416 outliers_IDs_list = []
2417 2417
2418 2418
2419 2419 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
2420 2420 'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels','cohFactor',
2421 2421 '__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'debug','prev_pnoise','thfactor')
2422 2422 def __init__(self, **kwargs):
2423 2423
2424 2424 Operation.__init__(self, **kwargs)
2425 2425 self.isConfig = False
2426 2426 self.currentTime = None
2427 2427
2428 2428 def setup(self,dataOut, n=None , navg=0.9, profileMargin=50,thHistOutlier=15,minHei=None, maxHei=None, nBins=10,
2429 2429 minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
2430 2430 idate=None,startH=None,endH=None, thfactor=1 ):
2431 2431
2432 2432 if n == None and timeInterval == None:
2433 2433 raise ValueError("nprofiles or timeInterval should be specified ...")
2434 2434
2435 2435 if n != None:
2436 2436 self.n = n
2437 2437
2438 2438 self.navg = navg
2439 2439 self.profileMargin = profileMargin
2440 2440 self.thHistOutlier = thHistOutlier
2441 2441 self.__profIndex = 0
2442 2442 self.buffer = None
2443 2443 self._ipp = dataOut.ippSeconds
2444 2444 self.n_prof_released = 0
2445 2445 self.heightList = dataOut.heightList
2446 2446 self.init_prof = 0
2447 2447 self.end_prof = 0
2448 2448 self.__count_exec = 0
2449 2449 self.__profIndex = 0
2450 2450 self.first_utcBlock = None
2451 2451 self.prev_pnoise = None
2452 2452 self.nBins = nBins
2453 2453 self.thfactor = thfactor
2454 2454 #self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
2455 2455 minHei = minHei
2456 2456 maxHei = maxHei
2457 2457 if minHei==None :
2458 2458 minHei = dataOut.heightList[0]
2459 2459 if maxHei==None :
2460 2460 maxHei = dataOut.heightList[-1]
2461 2461 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
2462 2462 self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
2463 2463 self.nChannels = dataOut.nChannels
2464 2464 self.nHeights = dataOut.nHeights
2465 2465 self.test_counter = 0
2466 2466 self.debug = debug
2467 2467 self.remYagi = remYagi
2468 2468 self.cohFactor = dataOut.nCohInt
2469 2469 if self.remYagi :
2470 2470 if minHJULIA==None or maxHJULIA==None:
2471 2471 raise ValueError("Parameters minHYagi and minHYagi are necessary!")
2472 2472 return
2473 2473 if idate==None or startH==None or endH==None:
2474 2474 raise ValueError("Date and hour parameters are necessary!")
2475 2475 return
2476 2476 self.minHJULIA_idx,self.maxHJULIA_idx = getHei_index(minHJULIA, maxHJULIA, dataOut.heightList)
2477 2477 self.offYagi = offYagi
2478 2478 self.nTxYagi = nProfYagi
2479 2479
2480 2480 self.startTime = datetime.datetime.combine(idate,startH)
2481 2481 self.endTime = datetime.datetime.combine(idate,endH)
2482 2482
2483 2483 log.warning("Be careful with the selection of parameters for sats removal! It is avisable to \
2484 2484 activate the debug parameter in this operation for calibration", self.name)
2485 2485
2486 2486
2487 2487 def filterSatsProfiles(self):
2488 2488
2489 2489 data = self.__buffer_data.copy()
2490 2490 #print(data.shape)
2491 2491 nChannels, profiles, heights = data.shape
2492 2492 indexes=numpy.zeros([], dtype=int)
2493 2493 indexes = numpy.delete(indexes,0)
2494 2494
2495 2495 indexesYagi=numpy.zeros([], dtype=int)
2496 2496 indexesYagi = numpy.delete(indexesYagi,0)
2497 2497
2498 2498 indexesYagi_up=numpy.zeros([], dtype=int)
2499 2499 indexesYagi_up = numpy.delete(indexesYagi_up,0)
2500 2500 indexesYagi_down=numpy.zeros([], dtype=int)
2501 2501 indexesYagi_down = numpy.delete(indexesYagi_down,0)
2502 2502
2503 2503
2504 2504 indexesJULIA=numpy.zeros([], dtype=int)
2505 2505 indexesJULIA = numpy.delete(indexesJULIA,0)
2506 2506
2507 2507 outliers_IDs=[]
2508 2508
2509 2509 div = profiles//self.nBins
2510 2510
2511 2511 for c in range(nChannels):
2512 2512 #print(self.min_ref,self.max_ref)
2513 2513
2514 2514 import scipy.signal
2515 2515 b, a = scipy.signal.butter(3, 0.5)
2516 2516 #noise_ref = (data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref]))
2517 2517 noise_ref = numpy.abs(data[c,:,self.min_ref:self.max_ref])
2518 2518 lnoise = len(noise_ref[0,:])
2519 2519 #print(noise_ref.shape)
2520 2520 noise_ref = noise_ref.mean(axis=1)
2521 2521 #fnoise = noise_ref
2522 2522 fnoise = scipy.signal.filtfilt(b, a, noise_ref)
2523 2523 #noise_refdB = 10* numpy.log10(noise_ref)
2524 2524 #print("Noise ",numpy.percentile(noise_ref,95))
2525 2525 p95 = numpy.percentile(fnoise,95)
2526 2526 mean_noise = fnoise.mean()
2527 2527
2528 2528 if self.prev_pnoise != None:
2529 2529 if mean_noise < (1.1 * self.prev_pnoise) and mean_noise > (0.9 * self.prev_pnoise):
2530 2530 mean_noise = 0.9*mean_noise + 0.1*self.prev_pnoise
2531 2531 self.prev_pnoise = mean_noise
2532 2532 else:
2533 2533 mean_noise = self.prev_pnoise
2534 2534 else:
2535 2535 self.prev_pnoise = mean_noise
2536 2536
2537 2537 std = fnoise.std()+ fnoise.mean()
2538 2538
2539 2539
2540 2540
2541 2541 #power = (data[c,:,self.minHei_idx:self.maxHei_idx] * numpy.conjugate(data[c,:,self.minHei_idx:self.maxHei_idx]))
2542 2542 power = numpy.abs(data[c,:,self.minHei_idx:self.maxHei_idx])
2543 2543 npower = len(power[0,:])
2544 2544 #print(power.shape)
2545 2545 power = power.mean(axis=1)
2546 2546
2547 2547 fpower = scipy.signal.filtfilt(b, a, power)
2548 2548 #print(power.shape)
2549 2549 #powerdB = 10* numpy.log10(power)
2550 2550
2551 2551 #th = p95 * self.thfactor
2552 2552 th = mean_noise * self.thfactor
2553 2553
2554 2554 index = numpy.where(fpower > th )
2555 2555 #print("Noise ",mean_noise, p95)
2556 2556 #print(index)
2557 2557
2558 2558
2559 2559
2560 2560 if index[0].size <= int(self.navg*profiles): #outliers from sats
2561 2561 indexes = numpy.append(indexes, index[0])
2562 2562
2563 2563 index2low = numpy.where(fpower < (th*0.5 )) #outliers from no TX
2564 2564 if index2low[0].size <= int(self.navg*profiles):
2565 2565 indexes = numpy.append(indexes, index2low[0])
2566 2566
2567 2567 #print("sdas ", noise_ref.mean())
2568 2568
2569 2569 if self.remYagi :
2570 2570 #print(self.minHJULIA_idx, self.maxHJULIA_idx)
2571 2571 powerJULIA = (data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx] * numpy.conjugate(data[c,:,self.minHJULIA_idx:self.maxHJULIA_idx])).real
2572 2572 powerJULIA = powerJULIA.mean(axis=1)
2573 2573 th_JULIA = powerJULIA.mean()*0.85
2574 2574 indexJULIA = numpy.where(powerJULIA >= th_JULIA )
2575 2575
2576 2576 indexesJULIA= numpy.append(indexesJULIA, indexJULIA[0])
2577 2577
2578 2578 # fig, ax = plt.subplots()
2579 2579 # ax.plot(powerJULIA)
2580 2580 # ax.axhline(th_JULIA, color='r')
2581 2581 # plt.grid()
2582 2582 # plt.show()
2583 2583
2584 2584 if self.debug:
2585 2585 fig, ax = plt.subplots()
2586 2586 ax.plot(fpower, label="power")
2587 2587 #ax.plot(fnoise, label="noise ref")
2588 2588 ax.axhline(th, color='g', label="th")
2589 2589 #ax.axhline(std, color='b', label="mean")
2590 2590 ax.legend()
2591 2591 plt.grid()
2592 2592 plt.show()
2593 2593
2594 2594 #print(indexes)
2595 2595
2596 2596 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
2597 2597 #outliers_IDs = numpy.unique(outliers_IDs)
2598 2598 # print(indexesJULIA)
2599 2599 if len(indexesJULIA > 1):
2600 2600 iJ = indexesJULIA
2601 2601 locs = [ (iJ[n]-iJ[n-1]) > 5 for n in range(len(iJ))]
2602 2602 locs_2 = numpy.where(locs)[0]
2603 2603 #print(locs_2, indexesJULIA[locs_2-1])
2604 2604 indexesYagi_up = numpy.append(indexesYagi_up, indexesJULIA[locs_2-1])
2605 2605 indexesYagi_down = numpy.append(indexesYagi_down, indexesJULIA[locs_2])
2606 2606
2607 2607
2608 2608 indexesYagi_up = numpy.append(indexesYagi_up,indexesJULIA[-1])
2609 2609 indexesYagi_down = numpy.append(indexesYagi_down,indexesJULIA[0])
2610 2610
2611 2611 indexesYagi_up = numpy.unique(indexesYagi_up)
2612 2612 indexesYagi_down = numpy.unique(indexesYagi_down)
2613 2613
2614 2614
2615 2615 aux_ind = [ numpy.arange( (self.offYagi + k)+1, (self.offYagi + k + self.nTxYagi)+1, 1, dtype=int) for k in indexesYagi_up]
2616 2616 indexesYagi_up = (numpy.asarray(aux_ind)).flatten()
2617 2617
2618 2618 aux_ind2 = [ numpy.arange( (k - self.nTxYagi)+1, k+1 , 1, dtype=int) for k in indexesYagi_down]
2619 2619 indexesYagi_down = (numpy.asarray(aux_ind2)).flatten()
2620 2620
2621 2621 indexesYagi = numpy.append(indexesYagi,indexesYagi_up)
2622 2622 indexesYagi = numpy.append(indexesYagi,indexesYagi_down)
2623 2623
2624 2624
2625 2625 indexesYagi = indexesYagi[ (indexesYagi >= 0) & (indexesYagi<profiles)]
2626 2626 indexesYagi = numpy.unique(indexesYagi)
2627 2627
2628 2628 #print("indexes: " ,indexes)
2629 2629 outs_lines = numpy.unique(indexes)
2630 2630 #print(outs_lines)
2631 2631
2632 2632 #Agrupando el histograma de outliers,
2633 2633 my_bins = numpy.linspace(0,int(profiles), div, endpoint=True)
2634 2634
2635 2635
2636 2636 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
2637 2637 #print("hist: ",hist)
2638 2638 hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier)[0] #es outlier
2639 2639 # print(hist_outliers_indexes)
2640 2640 if len(hist_outliers_indexes>0):
2641 2641 hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
2642 2642
2643 2643 bins_outliers_indexes = [int(i)+1 for i in (bins[hist_outliers_indexes])] #
2644 2644 outlier_loc_index = []
2645 2645 #print("out indexes ", bins_outliers_indexes)
2646 2646
2647 2647 # if len(bins_outliers_indexes) <= 2:
2648 2648 # extprof = 0
2649 2649 # else:
2650 2650 # extprof = self.profileMargin
2651 2651
2652 2652 extprof = self.profileMargin
2653 2653
2654 2654 outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-extprof,bins_outliers_indexes[n] + extprof) ]
2655 2655 outlier_loc_index = numpy.asarray(outlier_loc_index)
2656 2656 # if len(outlier_loc_index)>1:
2657 2657 # ipmax = numpy.where(fpower==fpower.max())[0]
2658 2658 # print("pmax: ",ipmax)
2659 2659
2660 2660
2661 2661
2662 2662
2663 2663 #print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
2664 2664 outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
2665 2665 #print("outliers final: ", outlier_loc_index)
2666 2666
2667 2667
2668 2668 if self.debug:
2669 2669 x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
2670 2670 fig, ax = plt.subplots(nChannels,2,figsize=(8, 6))
2671 2671
2672 2672 for i in range(nChannels):
2673 2673 dat = data[i,:,:].real
2674 2674 dat = 10* numpy.log10((data[i,:,:] * numpy.conjugate(data[i,:,:])).real)
2675 2675 m = numpy.nanmean(dat)
2676 2676 o = numpy.nanstd(dat)
2677 2677 if nChannels>1:
2678 2678 c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
2679 2679 ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
2680 2680 #fig.colorbar(c)
2681 2681 ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
2682 2682 ax[i][1].hist(outs_lines,bins=my_bins)
2683 2683 if self.remYagi :
2684 2684 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
2685 2685 else:
2686 2686 c = ax[0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = (70+2*self.cohFactor))
2687 2687 ax[0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
2688 2688 #fig.colorbar(c)
2689 2689 ax[0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
2690 2690
2691 2691 ax[1].hist(outs_lines,bins=my_bins)
2692 2692 if self.remYagi :
2693 2693 ax[0].vlines(indexesYagi,750,850, linestyles='dashed', label = 'yagi', color='m')
2694 2694 plt.show()
2695 2695
2696 2696
2697 2697
2698 2698
2699 2699 if self.remYagi and (self.currentTime < self.startTime and self.currentTime < self.endTime):
2700 2700 outlier_loc_index = numpy.append(outlier_loc_index,indexesYagi)
2701 2701
2702 2702 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2703 2703
2704 2704 #print("outs list: ", self.outliers_IDs_list)
2705 2705 return self.__buffer_data
2706 2706
2707 2707
2708 2708
2709 2709 def fillBuffer(self, data, datatime):
2710 2710
2711 2711 if self.__profIndex == 0:
2712 2712 self.__buffer_data = data.copy()
2713 2713
2714 2714 else:
2715 2715 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2716 2716 self.__profIndex += 1
2717 2717 self.__buffer_times.append(datatime)
2718 2718
2719 2719 def getData(self, data, datatime=None):
2720 2720
2721 2721 if self.__profIndex == 0:
2722 2722 self.__initime = datatime
2723 2723
2724 2724
2725 2725 self.__dataReady = False
2726 2726
2727 2727 self.fillBuffer(data, datatime)
2728 2728 dataBlock = None
2729 2729
2730 2730 if self.__profIndex == self.n:
2731 2731 #print("apnd : ",data)
2732 2732 dataBlock = self.filterSatsProfiles()
2733 2733 self.__dataReady = True
2734 2734
2735 2735 return dataBlock
2736 2736
2737 2737 if dataBlock is None:
2738 2738 return None, None
2739 2739
2740 2740
2741 2741
2742 2742 return dataBlock
2743 2743
2744 2744 def releaseBlock(self):
2745 2745
2746 2746 if self.n % self.lenProfileOut != 0:
2747 2747 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2748 2748 return None
2749 2749
2750 2750 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2751 2751
2752 2752 self.init_prof = self.end_prof
2753 2753 self.end_prof += self.lenProfileOut
2754 2754 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2755 2755 self.n_prof_released += 1
2756 2756
2757 2757 return data
2758 2758
2759 2759 def run(self, dataOut, n=None, navg=0.9, nProfilesOut=1, profile_margin=50, th_hist_outlier=15,minHei=None,nBins=10,
2760 2760 maxHei=None, minRef=None, maxRef=None, debug=False, remYagi=False, nProfYagi = 0, offYagi=0, minHJULIA=None, maxHJULIA=None,
2761 2761 idate=None,startH=None,endH=None, thfactor=1):
2762 2762
2763 2763 if not self.isConfig:
2764 2764 #print("init p idx: ", dataOut.profileIndex )
2765 2765 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,minHei=minHei,
2766 2766 nBins=10, maxHei=maxHei, minRef=minRef, maxRef=maxRef, debug=debug, remYagi=remYagi, nProfYagi = nProfYagi,
2767 2767 offYagi=offYagi, minHJULIA=minHJULIA,maxHJULIA=maxHJULIA,idate=idate,startH=startH,endH=endH, thfactor=thfactor)
2768 2768
2769 2769 self.isConfig = True
2770 2770
2771 2771 dataBlock = None
2772 2772 self.currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
2773 2773
2774 2774 if not dataOut.buffer_empty: #hay datos acumulados
2775 2775
2776 2776 if self.init_prof == 0:
2777 2777 self.n_prof_released = 0
2778 2778 self.lenProfileOut = nProfilesOut
2779 2779 dataOut.flagNoData = False
2780 2780 #print("tp 2 ",dataOut.data.shape)
2781 2781
2782 2782 self.init_prof = 0
2783 2783 self.end_prof = self.lenProfileOut
2784 2784
2785 2785 dataOut.nProfiles = self.lenProfileOut
2786 2786 if nProfilesOut == 1:
2787 2787 dataOut.flagDataAsBlock = False
2788 2788 else:
2789 2789 dataOut.flagDataAsBlock = True
2790 2790 #print("prof: ",self.init_prof)
2791 2791 dataOut.flagNoData = False
2792 2792 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
2793 2793 #print("omitting: ", self.n_prof_released)
2794 2794 dataOut.flagNoData = True
2795 2795 dataOut.ippSeconds = self._ipp
2796 2796 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
2797 2797 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
2798 2798 #dataOut.data = self.releaseBlock()
2799 2799 #########################################################3
2800 2800 if self.n % self.lenProfileOut != 0:
2801 2801 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2802 2802 return None
2803 2803
2804 2804 dataOut.data = None
2805 2805
2806 2806 if nProfilesOut == 1:
2807 2807 dataOut.data = self.buffer[:,self.end_prof-1,:] #ch, prof, alt
2808 2808 else:
2809 2809 dataOut.data = self.buffer[:,self.init_prof:self.end_prof,:] #ch, prof, alt
2810 2810
2811 2811 self.init_prof = self.end_prof
2812 2812 self.end_prof += self.lenProfileOut
2813 2813 #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
2814 2814 self.n_prof_released += 1
2815 2815
2816 2816 if self.end_prof >= (self.n +self.lenProfileOut):
2817 2817
2818 2818 self.init_prof = 0
2819 2819 self.__profIndex = 0
2820 2820 self.buffer = None
2821 2821 dataOut.buffer_empty = True
2822 2822 self.outliers_IDs_list = []
2823 2823 self.n_prof_released = 0
2824 2824 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
2825 2825 #print("cleaning...", dataOut.buffer_empty)
2826 2826 dataOut.profileIndex = self.__profIndex
2827 2827 ####################################################################
2828 2828 return dataOut
2829 2829
2830 2830
2831 2831 #print("tp 223 ",dataOut.data.shape)
2832 2832 dataOut.flagNoData = True
2833 2833
2834 2834
2835 2835
2836 2836 try:
2837 2837 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
2838 2838 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
2839 2839 self.__count_exec +=1
2840 2840 except Exception as e:
2841 2841 print("Error getting profiles data",self.__count_exec )
2842 2842 print(e)
2843 2843 sys.exit()
2844 2844
2845 2845 if self.__dataReady:
2846 2846 #print("omitting: ", len(self.outliers_IDs_list))
2847 2847 self.__count_exec = 0
2848 2848 #dataOut.data =
2849 2849 #self.buffer = numpy.flip(dataBlock, axis=1)
2850 2850 self.buffer = dataBlock
2851 2851 self.first_utcBlock = self.__initime
2852 2852 dataOut.utctime = self.__initime
2853 2853 dataOut.nProfiles = self.__profIndex
2854 2854 #dataOut.flagNoData = False
2855 2855 self.init_prof = 0
2856 2856 self.__profIndex = 0
2857 2857 self.__initime = None
2858 2858 dataBlock = None
2859 2859 self.__buffer_times = []
2860 2860 dataOut.error = False
2861 2861 dataOut.useInputBuffer = True
2862 2862 dataOut.buffer_empty = False
2863 2863 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
2864 2864
2865 2865
2866 2866
2867 2867 #print(self.__count_exec)
2868 2868
2869 2869 return dataOut
2870 2870
2871 2871
2872 2872
2873 2873
2874 2874 class remHeightsIppInterf(Operation):
2875 2875
2876 2876 def __init__(self, **kwargs):
2877 2877
2878 2878
2879 2879 Operation.__init__(self, **kwargs)
2880 2880
2881 2881 self.isConfig = False
2882 2882
2883 2883 self.heights_indx = None
2884 2884 self.heightsList = []
2885 2885
2886 2886 self.ipp1 = None
2887 2887 self.ipp2 = None
2888 2888 self.tx1 = None
2889 2889 self.tx2 = None
2890 2890 self.dh1 = None
2891 2891
2892 2892
2893 2893 def setup(self, dataOut, ipp1=None, ipp2=None, tx1=None, tx2=None, dh1=None,
2894 2894 idate=None, startH=None, endH=None):
2895 2895
2896 2896
2897 2897 self.ipp1 = ipp1
2898 2898 self.ipp2 = ipp2
2899 2899 self.tx1 = tx1
2900 2900 self.tx2 = tx2
2901 2901 self.dh1 = dh1
2902 2902
2903 2903 _maxIpp1R = dataOut.heightList.max()
2904 2904
2905 2905 _n_repeats = int(_maxIpp1R / ipp2)
2906 2906 _init_hIntf = (tx1 + ipp2/2)+ dh1
2907 2907 _n_hIntf = int(tx2 / dh1)
2908 2908
2909 2909 self.heightsList = [_init_hIntf+n*ipp2 for n in range(_n_repeats) ]
2910 2910 heiList = dataOut.heightList
2911 2911 self.heights_indx = [getHei_index(h,h,heiList)[0] for h in self.heightsList]
2912 2912
2913 2913 self.heights_indx = [ numpy.asarray([k for k in range(_n_hIntf+2)])+(getHei_index(h,h,heiList)[0] -1) for h in self.heightsList]
2914 2914
2915 2915 self.heights_indx = numpy.asarray(self.heights_indx )
2916 2916 self.isConfig = True
2917 2917 self.startTime = datetime.datetime.combine(idate,startH)
2918 2918 self.endTime = datetime.datetime.combine(idate,endH)
2919 2919 #print(self.startTime, self.endTime)
2920 2920 #print("nrepeats: ", _n_repeats, " _nH: ",_n_hIntf )
2921 2921
2922 2922 log.warning("Heights set to zero (km): ", self.name)
2923 2923 log.warning(str((dataOut.heightList[self.heights_indx].flatten())), self.name)
2924 2924 log.warning("Be careful with the selection of heights for noise calculation!")
2925 2925
2926 2926 def run(self, dataOut, ipp1=None, ipp2=None, tx1=None, tx2=None, dh1=None, idate=None,
2927 2927 startH=None, endH=None):
2928 2928 #print(locals().values())
2929 2929 if None in locals().values():
2930 2930 log.warning('Missing kwargs, invalid values """None""" ', self.name)
2931 2931 return dataOut
2932 2932
2933 2933
2934 2934 if not self.isConfig:
2935 2935 self.setup(dataOut, ipp1=ipp1, ipp2=ipp2, tx1=tx1, tx2=tx2, dh1=dh1,
2936 2936 idate=idate, startH=startH, endH=endH)
2937 2937
2938 2938 dataOut.flagProfilesByRange = False
2939 2939 currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
2940 2940
2941 2941 if currentTime < self.startTime or currentTime > self.endTime:
2942 2942 return dataOut
2943 2943
2944 2944 for ch in range(dataOut.data.shape[0]):
2945 2945
2946 2946 for hk in self.heights_indx.flatten():
2947 2947 if dataOut.data.ndim < 3:
2948 2948 dataOut.data[ch,hk] = 0.0 + 0.0j
2949 2949 else:
2950 2950 dataOut.data[ch,:,hk] = 0.0 + 0.0j
2951 2951
2952 2952 dataOut.flagProfilesByRange = True
2953 2953
2954 2954 return dataOut
2955 2955
2956 2956
2957 2957
2958 2958
2959 2959 class profiles2Block(Operation):
2960 2960 '''
2961 2961 Escrito: Joab Apaza
2962 2962
2963 2963 genera un bloque de perfiles
2964 2964
2965 2965
2966 2966 Out:
2967 2967 block
2968 2968 '''
2969 2969
2970 2970 isConfig = False
2971 2971 __buffer_data = []
2972 2972 __buffer_times = []
2973 2973 __profIndex = 0
2974 2974 __byTime = False
2975 2975 __initime = None
2976 2976 __lastdatatime = None
2977 2977 buffer = None
2978 2978 n = None
2979 2979 __dataReady = False
2980 2980 __nChannels = None
2981 2981 __nHeis = None
2982 2982
2983 2983 def __init__(self, **kwargs):
2984 2984
2985 2985 Operation.__init__(self, **kwargs)
2986 2986 self.isConfig = False
2987 2987
2988 2988 def setup(self,n=None, timeInterval=None):
2989 2989
2990 2990 if n == None and timeInterval == None:
2991 2991 raise ValueError("n or timeInterval should be specified ...")
2992 2992
2993 2993 if n != None:
2994 2994 self.n = n
2995 2995 self.__byTime = False
2996 2996 else:
2997 2997 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
2998 2998 self.n = 9999
2999 2999 self.__byTime = True
3000 3000
3001 3001 self.__profIndex = 0
3002 3002
3003 3003
3004 3004 def fillBuffer(self, data, datatime):
3005 3005
3006 3006 if self.__profIndex == 0:
3007 3007 self.__buffer_data = data.copy()
3008 3008
3009 3009 else:
3010 3010 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
3011 3011 self.__profIndex += 1
3012 3012 self.__buffer_times.append(datatime)
3013 3013
3014 3014 def getData(self, data, datatime=None):
3015 3015 if self.__initime == None:
3016 3016 self.__initime = datatime
3017 3017
3018 3018 if data.ndim < 3:
3019 3019 data = data.reshape(self.__nChannels,1,self.__nHeis )
3020 3020
3021 3021 if self.__byTime:
3022 3022 dataBlock = self.byTime(data, datatime)
3023 3023 else:
3024 3024 dataBlock = self.byProfiles(data, datatime)
3025 3025
3026 3026
3027 3027 self.__lastdatatime = datatime
3028 3028
3029 3029 if dataBlock is None:
3030 3030 return None, None
3031 3031
3032 3032 return dataBlock, self.__buffer_times
3033 3033
3034 3034 def byProfiles(self, data, datatime):
3035 3035
3036 3036 self.__dataReady = False
3037 3037 dataBlock = None
3038 3038 # n = None
3039 3039 # print data
3040 3040 # raise
3041 3041 self.fillBuffer(data, datatime)
3042 3042
3043 3043 if self.__profIndex == self.n:
3044 3044 dataBlock = self.__buffer_data
3045 3045 self.__dataReady = True
3046 3046
3047 3047 return dataBlock
3048 3048
3049 3049 def byTime(self, data, datatime):
3050 3050
3051 3051 self.__dataReady = False
3052 3052 dataBlock = None
3053 3053 n = None
3054 3054
3055 3055 self.fillBuffer(data, datatime)
3056 3056
3057 3057 if (datatime - self.__initime) >= self.__integrationtime:
3058 3058 dataBlock = self.__buffer_data
3059 3059 self.n = self.__profIndex
3060 3060 self.__dataReady = True
3061 3061
3062 3062 return dataBlock
3063 3063
3064 3064
3065 3065 def run(self, dataOut, n=None, timeInterval=None, **kwargs):
3066 3066
3067 3067 if not self.isConfig:
3068 3068 self.setup(n=n, timeInterval=timeInterval, **kwargs)
3069 3069 self.__nChannels = dataOut.nChannels
3070 3070 self.__nHeis = len(dataOut.heightList)
3071 3071 self.isConfig = True
3072 3072
3073 3073 if dataOut.flagDataAsBlock:
3074 3074 """
3075 3075 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
3076 3076 """
3077 3077 raise ValueError("The data is already a block")
3078 3078 return
3079 3079 else:
3080 3080
3081 3081 dataBlock, timeBlock = self.getData(dataOut.data, dataOut.utctime)
3082 3082
3083 3083
3084 3084 # print(dataOut.data.shape)
3085 3085 # dataOut.timeInterval *= n
3086 3086 dataOut.flagNoData = True
3087 3087
3088 3088 if self.__dataReady:
3089 3089 dataOut.data = dataBlock
3090 3090 dataOut.flagDataAsBlock = True
3091 3091 dataOut.utctime = timeBlock[-1]
3092 3092 dataOut.nProfiles = self.__profIndex
3093 3093 # print avgdata, avgdatatime
3094 3094 # raise
3095 3095 dataOut.flagNoData = False
3096 3096 self.__profIndex = 0
3097 3097 self.__initime = None
3098 3098 #update Processing Header:
3099 3099 # print(dataOut.data.shape)
3100 3100
3101 3101 return dataOut
3102 3102
3103 3103
3104 3104 class remFaradayProfiles(Operation):
3105 3105
3106 3106 def __init__(self, **kwargs):
3107 3107
3108 3108
3109 3109 Operation.__init__(self, **kwargs)
3110 3110
3111 3111 self.isConfig = False
3112 3112
3113 3113 self.nprofile2 = 0
3114 3114 self.profile = 0
3115 3115 self.flagRun = False
3116 3116 self.flagRemove = False
3117 3117 self.k = 0
3118 3118
3119 def setup(self, channel,nChannels=5, nProfiles=300,nBlocks=100, nIpp2=200, nTx2=132, nTaus=22, offTaus=14, iTaus=8,
3119 def setup(self, channel,nChannels=5, nProfiles=300,nBlocks=100, nIpp2=300, nTx2=132, nTaus=22, offTaus=14, iTaus=8,
3120 3120 nfft=1):
3121 3121 '''
3122 3122 nProfiles = amisr profiles per block -> raw data
3123 3123 nIpp1 = number of profiles in one AMISR sync
3124 3124 nIpp2 = number of profiles in one Jicamarca sync
3125 3125 nTx2 = number of profiles transmited for Faraday Experiment
3126 3126 nTaus = Total profiles for lags
3127 3127 offTaus = where starts the interference, (profile)
3128 3128 iTaus = lenght of the interference
3129 3129 irepeat = number of repetition of the Taus
3130 3130 '''
3131 3131 self.nIpp2 = nIpp2
3132 3132 self.channel = channel
3133 3133 self.nChannels = nChannels
3134 3134 self.nTx2 = nTx2
3135 3135 self.nTaus = nTaus
3136 3136
3137 3137
3138 3138 booldataset = numpy.ones( (nBlocks, nProfiles) )
3139 3139 self.profilesFlag = None
3140 3140 #marking the afected profiles
3141 3141 f_iTaus=False
3142 3142 f_ntx = False
3143 3143 fi = 0
3144 3144 k = 0
3145 3145 kt =0
3146 3146 fi_reps = 0
3147 3147 for i in range(nBlocks):
3148 3148 for j in range(nProfiles):
3149 3149 # fi 0---nTaus
3150 3150 #
3151 3151 if k%nIpp2==0: #each sync PPs or 2, 3, or 5
3152 3152 f_ntx = True
3153 3153 kt = 0
3154 #print(k, fi, j, f_iTaus)
3155 3154 if f_ntx:
3156 3155
3157 3156 if kt%nTaus==0: #each sequence of Taus
3158 3157 f_iTaus = True
3159 3158 fi = 0
3160 3159
3161 3160 if f_iTaus:
3162 3161 if fi > offTaus-1:
3163 3162 booldataset[i, j]=0 #Afected profile
3164 3163 fi += 1
3165 3164 if fi == nTaus-1: #restart the taus sequence
3166 3165 fi = 0
3167 3166 f_iTaus = False
3168 3167 fi_reps += 1
3169 3168 if fi_reps == (nTx2/nTaus):
3170 3169 fi = 0
3171 #print("AQUI, ", fi_reps, k, fi)
3172 3170 fi_reps = 0
3173 3171 f_ntx=False
3174 # if i < 1:
3175 # print(fi, kt)
3176 3172 kt += 1
3177 3173 k += 1
3178 3174
3179 3175 # fig = plt.figure()
3180 3176 # ax = fig.add_subplot(111)
3181 3177 # cax = ax.pcolormesh(booldataset, cmap='plasma')
3182 3178 # cbar = fig.colorbar(cax)
3183 3179 # plt.show()
3184 3180
3185
3186 #print ("AQUI")
3181
3187 3182 #reshape the Flag as AMISR reader
3188 3183
3189 3184 profPerCH = int( (nProfiles) / (nfft*nChannels))
3190 3185 new_block = numpy.empty( (nBlocks, nChannels, int(nProfiles/nChannels) ) )
3191 3186 # print(new_block.shape, profPerCH)
3192 3187 for thisChannel in range(nChannels):
3193 3188
3194 3189 ich = thisChannel
3195 3190
3196 3191 idx_ch = [nfft*(ich + nChannels*k) for k in range(profPerCH)]
3197 3192 #print(idx_ch)
3198 3193 if nfft > 1:
3199 3194 aux = [numpy.arange(i, i+nfft) for i in idx_ch]
3200 3195 idx_ch = None
3201 3196 idx_ch =aux
3202 3197 idx_ch = numpy.array(idx_ch, dtype=int).flatten()
3203 3198 else:
3204 3199 idx_ch = numpy.array(idx_ch, dtype=int)
3205 3200
3206 3201 new_block[:,ich,:] = booldataset[:,idx_ch]
3207 3202
3208 3203 new_block = numpy.transpose(new_block, (1,0,2))
3209 new_block = numpy.reshape(new_block, (nChannels,-1))
3210 #new_block = numpy.reshape(new_block, (nChannels,profPerCH*nBlocks))
3204 #new_block = numpy.reshape(new_block, (nChannels,-1))
3205 new_block = numpy.reshape(new_block, (nChannels,profPerCH*nBlocks))
3211 3206 self.profilesFlag = new_block.copy()
3212 3207
3213 3208 # fig = plt.figure()
3214 3209 # ax = fig.add_subplot(111)
3215 3210 # cax = ax.pcolormesh(new_block, cmap='plasma')
3216 3211 # cbar = fig.colorbar(cax)
3217 3212 # plt.show()
3218 3213
3219 3214 self.isConfig = True
3220 3215
3221 #print(self.profilesFlag.shape)
3222 3216
3223 3217 def run(self,dataOut, channel=0, nChannels=5, nProfiles=300,nBlocks=100,nIpp1=100,
3224 3218 nIpp2=300, nTx2=132, nTaus=22, offTaus=8, iTaus=14, nfft=1 ,offIpp=0):
3225 3219
3226 3220 dataOut.flagNoData = False
3227 3221
3228 3222 if not self.isConfig:
3229 3223 self.setup(channel,nChannels=nChannels, nProfiles=nProfiles,nBlocks=nBlocks, nIpp2=nIpp2,
3230 3224 nTx2=nTx2, nTaus=nTaus, offTaus=offTaus, iTaus=iTaus, nfft=nfft)
3231 3225 #print("Setup Done")
3232 3226 #print(offIpp*nIpp1/nChannels)
3233 3227 if not self.flagRun:
3234 3228 if self.nprofile2 < offIpp*nIpp1/nChannels :
3235 3229 self.nprofile2 += 1
3236 3230 return dataOut
3237 3231 else:
3238 3232 self.flagRun = True
3239 3233 self.profile = 0
3240
3241 3234
3242 3235 #check profile ## Faraday interference
3243 3236 if self.profilesFlag[channel, self.profile]==0:
3244 3237 dataOut.flagNoData = True # do not pass this profile
3245 # print(self.nprofile, dataOut.flagNoData)
3246 #print(self.nprofile2, self.profile, dataOut.flagNoData)
3238
3247 3239 self.profile +=1
3248 # if self.profile == int((nProfiles*nBlocks)/self.nChannels):
3249 # self.flagRun=False
3250 # self.profile = 0
3240
3251 3241
3252 3242 self.nprofile2 +=1
3253 3243
3254 3244 if self.nprofile2 == int((nProfiles*nBlocks)/self.nChannels):
3255 3245 self.nprofile2 = 0
3256 3246 self.profile = 0
3257 3247 self.flagRun = False
3258 3248
3259 3249 return dataOut No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now