##// END OF EJS Templates
schain
RSS
max_nIncohInt changed as @property, correction for normFactor used as matrix
joabAM -
r1554:7454ff2c4b3d
parent child
Show More
Side by Side Unified
Context file:
r1554:7454ff2c4b3d Loading diff...:
Show file before | Show file after | Hide comments
@@ -204,7 +204,7 class JROData(GenericData):
204 sampled_heightsFFT = None
204 sampled_heightsFFT = None
205 pulseLength_TxA = None
205 pulseLength_TxA = None
206 deltaHeight = None
206 deltaHeight = None
207
207
208 def __str__(self):
208 def __str__(self):
209
209
210 return '{} - {}'.format(self.type, self.datatime())
210 return '{} - {}'.format(self.type, self.datatime())
@@ -443,7 +443,7 class Voltage(JROData):
443
443
444 class Spectra(JROData):
444 class Spectra(JROData):
445
445
446 data_outlier = 0
446 data_outlier = None
447
447
448 def __init__(self):
448 def __init__(self):
449 '''
449 '''
@@ -483,7 +483,6 class Spectra(JROData):
483 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
483 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
484
484
485
485
486 self.max_nIncohInt = 1
487
486
488
487
489 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
488 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
@@ -511,6 +510,7 class Spectra(JROData):
511 #
510 #
512 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
511 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513 noise = numpy.zeros(self.nChannels)
512 noise = numpy.zeros(self.nChannels)
513
514 for channel in range(self.nChannels):
514 for channel in range(self.nChannels):
515 daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
515 daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
516
516
@@ -576,6 +576,7 class Spectra(JROData):
576
576
577 if self.flagDecodeData:
577 if self.flagDecodeData:
578 pwcode = numpy.sum(self.code[0]**2)
578 pwcode = numpy.sum(self.code[0]**2)
579 #print(self.flagDecodeData, pwcode)
579 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
580 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
580 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
581 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
581
582
@@ -612,15 +613,36 class Spectra(JROData):
612 factor = self.normFactor
613 factor = self.normFactor
613 power = numpy.zeros( (self.nChannels,self.nHeights) )
614 power = numpy.zeros( (self.nChannels,self.nHeights) )
614 for ch in range(self.nChannels):
615 for ch in range(self.nChannels):
616 z = None
615 if hasattr(factor,'shape'):
617 if hasattr(factor,'shape'):
616 z = numpy.divide(self.data_spc[ch],factor[ch])
618 if factor.ndim > 1:
619 z = self.data_spc[ch]/factor[ch]
620 else:
621 z = self.data_spc[ch]/factor
617 else:
622 else:
618 z = numpy.divide(self.data_spc[ch],factor)
623 z = self.data_spc[ch]/factor
619 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
624 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
620 avg = numpy.average(z, axis=0)
625 avg = numpy.average(z, axis=0)
621 power[ch,:] = 10 * numpy.log10(avg)
626 power[ch] = 10 * numpy.log10(avg)
622 return power
627 return power
623
628
629 @property
630 def max_nIncohInt(self):
631
632 ints = numpy.zeros(self.nChannels)
633 for ch in range(self.nChannels):
634 if hasattr(self.nIncohInt,'shape'):
635 if self.nIncohInt.ndim > 1:
636 ints[ch,] = self.nIncohInt[ch].max()
637 else:
638 ints[ch,] = self.nIncohInt
639 self.nIncohInt = int(self.nIncohInt)
640 else:
641 ints[ch,] = self.nIncohInt
642
643 return ints
644
645
624 def getCoherence(self, pairsList=None, phase=False):
646 def getCoherence(self, pairsList=None, phase=False):
625
647
626 z = []
648 z = []
@@ -901,7 +923,7 class Parameters(Spectra):
901 nAvg = None
923 nAvg = None
902 noise_estimation = None
924 noise_estimation = None
903 GauSPC = None # Fit gaussian SPC
925 GauSPC = None # Fit gaussian SPC
904 max_nIncohInt = 1
926
905 def __init__(self):
927 def __init__(self):
906 '''
928 '''
907 Constructor
929 Constructor
Show file before | Show file after | Hide comments
@@ -56,12 +56,25 class SpectraPlot(Plot):
56 data = {}
56 data = {}
57 meta = {}
57 meta = {}
58
58
59 data['rti'] = dataOut.getPower()
59 #data['rti'] = dataOut.getPower()
60 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
60 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
61 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
61 noise = 10*numpy.log10(dataOut.getNoise()/norm)
62 data['noise'] = noise
62
63 spc = 10*numpy.log10(dataOut.data_spc/norm)
63
64 z = []
65 for ch in range(dataOut.nChannels):
66 if hasattr(dataOut.normFactor,'shape'):
67 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
68 else:
69 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
70
71 z = numpy.asarray(z)
72 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
73 spc = 10*numpy.log10(z)
74
64 data['spc'] = spc
75 data['spc'] = spc
76 data['rti'] = spc.mean(axis=1)
77 data['noise'] = noise
65 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
78 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
66 if self.CODE == 'spc_moments':
79 if self.CODE == 'spc_moments':
67 data['moments'] = dataOut.moments
80 data['moments'] = dataOut.moments
@@ -270,10 +283,11 class RTIPlot(Plot):
270 self.update_list(dataOut)
283 self.update_list(dataOut)
271 data = {}
284 data = {}
272 meta = {}
285 meta = {}
286
273 data['rti'] = dataOut.getPower()
287 data['rti'] = dataOut.getPower()
274
288
275 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
289 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
276 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
290 noise = 10*numpy.log10(dataOut.getNoise()/norm)
277 data['noise'] = noise
291 data['noise'] = noise
278
292
279 return data, meta
293 return data, meta
@@ -529,12 +543,23 class SpectraCutPlot(Plot):
529 data = {}
543 data = {}
530 meta = {}
544 meta = {}
531
545
532 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
546 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
533 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
547 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
534
535 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
536 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
548 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
537
549
550
551 z = []
552 for ch in range(dataOut.nChannels):
553 if hasattr(dataOut.normFactor,'shape'):
554 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
555 else:
556 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
557
558 z = numpy.asarray(z)
559 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
560 spc = 10*numpy.log10(z)
561
562
538 data['spc'] = spc - noise
563 data['spc'] = spc - noise
539 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
564 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
540
565
@@ -860,18 +885,27 class NoiselessSpectraPlot(Plot):
860 data = {}
885 data = {}
861 meta = {}
886 meta = {}
862
887
863 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
888 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
864 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
889 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
890 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
865
891
892 z = []
893 for ch in range(dataOut.nChannels):
894 if hasattr(dataOut.normFactor,'shape'):
895 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
896 else:
897 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
866
898
867 #spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
899 z = numpy.asarray(z)
868 spc = 10*numpy.log10(dataOut.data_spc/norm)
900 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
901 spc = 10*numpy.log10(z)
869
902
870 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
871 data['rti'] = dataOut.getPower() - noise
872
903
873 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
874 data['spc'] = spc - noise
904 data['spc'] = spc - noise
905 #print(spc.shape)
906 data['rti'] = spc.mean(axis=1)
907 data['noise'] = noise
908
875
909
876
910
877 # data['noise'] = noise
911 # data['noise'] = noise
@@ -965,12 +999,12 class NoiselessRTIPlot(Plot):
965 self.update_list(dataOut)
999 self.update_list(dataOut)
966 data = {}
1000 data = {}
967 meta = {}
1001 meta = {}
1002 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1003 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt)
1004 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
968
1005
969
1006
970 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1007 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
971 #print("Norm: ", norm)
972 #print(dataOut.nProfiles , dataOut.max_nIncohInt ,dataOut.nCohInt, dataOut.windowOfFilter)
973 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
974
1008
975 data['noise'] = n0
1009 data['noise'] = n0
976 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1010 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
Show file before | Show file after | Hide comments
@@ -274,6 +274,7 class MergeH5(object):
274 elif isinstance(meta[name], dict):
274 elif isinstance(meta[name], dict):
275 for key, value in meta[name].items():
275 for key, value in meta[name].items():
276 return value[x]
276 return value[x]
277
277 if 'cspc' in name:
278 if 'cspc' in name:
278 return 'pair{:02d}'.format(x)
279 return 'pair{:02d}'.format(x)
279 else:
280 else:
@@ -315,12 +316,25 class MergeH5(object):
315 self.dataOut.utctime = time
316 self.dataOut.utctime = time
316 ints = [data.nIncohInt for data in self.ch_dataIn]
317 ints = [data.nIncohInt for data in self.ch_dataIn]
317 self.dataOut.nIncohInt = numpy.stack(ints, axis=1)
318 self.dataOut.nIncohInt = numpy.stack(ints, axis=1)
319
320 print("nIncohInt 1: ",self.dataOut.nIncohInt.shape)
321
322 if self.dataOut.nIncohInt.ndim > 3:
323 aux = self.dataOut.nIncohInt
324 self.dataOut.nIncohInt = None
325 self.dataOut.nIncohInt = aux[0]
326
327
318 if self.dataOut.nIncohInt.ndim < 3:
328 if self.dataOut.nIncohInt.ndim < 3:
319 nIncohInt = numpy.repeat(self.dataOut.nIncohInt, self.dataOut.nHeights).reshape(self.blocksPerFile,self.nChannels, self.dataOut.nHeights)
329 nIncohInt = numpy.repeat(self.dataOut.nIncohInt, self.dataOut.nHeights).reshape(self.blocksPerFile,self.nChannels, self.dataOut.nHeights)
320 #nIncohInt = numpy.reshape(nIncohInt, (self.blocksPerFile,self.nChannels, self.dataOut.nHeights))
330 #nIncohInt = numpy.reshape(nIncohInt, (self.blocksPerFile,self.nChannels, self.dataOut.nHeights))
321 self.dataOut.nIncohInt = None
331 self.dataOut.nIncohInt = None
322 self.dataOut.nIncohInt = nIncohInt
332 self.dataOut.nIncohInt = nIncohInt
323 #print("nIncohInt: ", self.dataOut.nIncohInt.shape)
333
334 if (self.dataOut.nIncohInt.shape)[0]==self.nChannels: ## ch,blocks, hei
335 self.dataOut.nIncohInt = numpy.swapaxes(self.dataOut.nIncohInt, 0, 1) ## blocks,ch, hei
336
337 print("nIncohInt 2: ", self.dataOut.nIncohInt.shape)
324 #print("utcTime: ", time.shape)
338 #print("utcTime: ", time.shape)
325 #print("data_spc ",self.dataOut.data_spc.shape)
339 #print("data_spc ",self.dataOut.data_spc.shape)
326 pairsList = [pair for pair in itertools.combinations(self.channelList, 2)]
340 pairsList = [pair for pair in itertools.combinations(self.channelList, 2)]
@@ -371,8 +385,10 class MergeH5(object):
371 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
385 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
372 dsDict['nDim'] = 0
386 dsDict['nDim'] = 0
373 else:
387 else:
388
374 dsDict['nDim'] = len(dataAux.shape) -1
389 dsDict['nDim'] = len(dataAux.shape) -1
375 dsDict['shape'] = dataAux.shape
390 dsDict['shape'] = dataAux.shape
391
376 if len(dsDict['shape'])>=2:
392 if len(dsDict['shape'])>=2:
377 dsDict['dsNumber'] = dataAux.shape[1]
393 dsDict['dsNumber'] = dataAux.shape[1]
378 else:
394 else:
@@ -399,6 +415,7 class MergeH5(object):
399 self.ch_dataIn[ch].utctime = None
415 self.ch_dataIn[ch].utctime = None
400 self.ch_dataIn[ch].nIncohInt = None
416 self.ch_dataIn[ch].nIncohInt = None
401 self.meta ={}
417 self.meta ={}
418 self.blocksPerFile = None
402
419
403 def writeData(self, outFilename):
420 def writeData(self, outFilename):
404
421
@@ -422,15 +439,17 class MergeH5(object):
422 sgrp = grp.create_group(label)
439 sgrp = grp.create_group(label)
423 else:
440 else:
424 sgrp = grp
441 sgrp = grp
442 k = -1*(dsInfo['nDim'] - 1)
443 #print(k, dsInfo['shape'], dsInfo['shape'][k:])
425 for i in range(dsInfo['dsNumber']):
444 for i in range(dsInfo['dsNumber']):
426 ds = sgrp.create_dataset(
445 ds = sgrp.create_dataset(
427 self.getLabel(dsInfo['variable'], i),
446 self.getLabel(dsInfo['variable'], i),(self.blocksPerFile, ) + dsInfo['shape'][k:],
428 (self.blocksPerFile, ) + dsInfo['shape'][2:],
429 chunks=True,
447 chunks=True,
430 dtype=dsInfo['dtype'])
448 dtype=dsInfo['dtype'])
431 dtsets.append(ds)
449 dtsets.append(ds)
432 data.append((dsInfo['variable'], i))
450 data.append((dsInfo['variable'], i))
433
451
452 #print("\n",dtsets)
434
453
435 print('Creating merged file: {}'.format(fp.filename))
454 print('Creating merged file: {}'.format(fp.filename))
436
455
@@ -442,6 +461,7 class MergeH5(object):
442 #print(ds, getattr(self.dataOut, attr)[ch].shape)
461 #print(ds, getattr(self.dataOut, attr)[ch].shape)
443 aux = getattr(self.dataOut, attr)# block, ch, ...
462 aux = getattr(self.dataOut, attr)# block, ch, ...
444 aux = numpy.swapaxes(aux,0,1) # ch, blocks, ...
463 aux = numpy.swapaxes(aux,0,1) # ch, blocks, ...
464 #print(ds.shape, aux.shape)
445 #ds[:] = getattr(self.dataOut, attr)[ch]
465 #ds[:] = getattr(self.dataOut, attr)[ch]
446 ds[:] = aux[ch]
466 ds[:] = aux[ch]
447
467
@@ -468,7 +488,7 class MergeH5(object):
468 self.readFile(fp,ch)
488 self.readFile(fp,ch)
469 fp.close()
489 fp.close()
470 self.getDataOut()
490 self.getDataOut()
471 name = name[-16:-1]
491 name = name[-16:]
472 #print("Final name out: ", name)
492 #print("Final name out: ", name)
473 outFile = os.path.join(self.pathOut, name)
493 outFile = os.path.join(self.pathOut, name)
474 #print("Outfile: ", outFile)
494 #print("Outfile: ", outFile)
Show file before | Show file after | Hide comments
@@ -141,7 +141,6 class ParametersProc(ProcessingUnit):
141 self.dataOut.nIncohInt = self.dataIn.nIncohInt
141 self.dataOut.nIncohInt = self.dataIn.nIncohInt
142 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
142 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
143 self.dataOut.ippFactor = self.dataIn.ippFactor
143 self.dataOut.ippFactor = self.dataIn.ippFactor
144 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
145 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
144 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
146 self.dataOut.radarControllerHeaderTxt = self.dataOut.radarControllerHeaderObj.toString()
145 self.dataOut.radarControllerHeaderTxt = self.dataOut.radarControllerHeaderObj.toString()
147 self.dataOut.ipp = self.dataIn.ipp
146 self.dataOut.ipp = self.dataIn.ipp
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (938 lines changed) Show them Hide them
@@ -63,7 +63,6 class SpectraProc(ProcessingUnit):
63 self.dataOut.flagShiftFFT = False
63 self.dataOut.flagShiftFFT = False
64 self.dataOut.nCohInt = self.dataIn.nCohInt
64 self.dataOut.nCohInt = self.dataIn.nCohInt
65 self.dataOut.nIncohInt = 1
65 self.dataOut.nIncohInt = 1
66 self.dataOut.max_nIncohInt = 1
67 self.dataOut.radar_ipp = self.dataIn.radar_ipp
66 self.dataOut.radar_ipp = self.dataIn.radar_ipp
68 self.dataOut.sampled_heightsFFT = self.dataIn.sampled_heightsFFT
67 self.dataOut.sampled_heightsFFT = self.dataIn.sampled_heightsFFT
69 self.dataOut.pulseLength_TxA = self.dataIn.pulseLength_TxA
68 self.dataOut.pulseLength_TxA = self.dataIn.pulseLength_TxA
@@ -690,11 +689,12 class getNoiseB(Operation):
690
689
691 self.dataOut.noise_estimation = None
690 self.dataOut.noise_estimation = None
692 noise = None
691 noise = None
692 #print("data type: ",self.dataOut.type, self.dataOut.nIncohInt, self.dataOut.max_nIncohInt)
693 if self.dataOut.type == 'Voltage':
693 if self.dataOut.type == 'Voltage':
694 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
694 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
695 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
695 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
696 elif self.dataOut.type == 'Spectra':
696 elif self.dataOut.type == 'Spectra':
697 #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt.shape)
697 #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt)
698 noise = numpy.zeros( self.dataOut.nChannels)
698 noise = numpy.zeros( self.dataOut.nChannels)
699 norm = 1
699 norm = 1
700
700
@@ -702,8 +702,8 class getNoiseB(Operation):
702 if not hasattr(self.dataOut.nIncohInt,'__len__'):
702 if not hasattr(self.dataOut.nIncohInt,'__len__'):
703 norm = 1
703 norm = 1
704 else:
704 else:
705 norm = self.dataOut.max_nIncohInt/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
705 norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
706 #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape)
706 #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape, self.dataOut.max_nIncohInt)
707 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
707 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
708 daux = numpy.multiply(daux, norm)
708 daux = numpy.multiply(daux, norm)
709 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
709 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
@@ -711,7 +711,8 class getNoiseB(Operation):
711 #print(daux.shape, daux)
711 #print(daux.shape, daux)
712 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
712 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
713 sortdata = numpy.sort(daux, axis=None)
713 sortdata = numpy.sort(daux, axis=None)
714 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset
714
715 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset
715
716
716
717
717 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
718 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
@@ -721,7 +722,7 class getNoiseB(Operation):
721 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
722 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
722 #print("2: ",self.dataOut.noise_estimation)
723 #print("2: ",self.dataOut.noise_estimation)
723 #print(self.dataOut.flagNoData)
724 #print(self.dataOut.flagNoData)
724 # print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor)
725 #print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor)
725 return self.dataOut
726 return self.dataOut
726
727
727 def getNoiseByMean(self,data):
728 def getNoiseByMean(self,data):
@@ -798,464 +799,464 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
798 return y
799 return y
799
800
800
801
801 class CleanRayleigh(Operation):
802 # class CleanRayleigh(Operation):
802
803 #
803 def __init__(self):
804 # def __init__(self):
804
805 #
805 Operation.__init__(self)
806 # Operation.__init__(self)
806 self.i=0
807 # self.i=0
807 self.isConfig = False
808 # self.isConfig = False
808 self.__dataReady = False
809 # self.__dataReady = False
809 self.__profIndex = 0
810 # self.__profIndex = 0
810 self.byTime = False
811 # self.byTime = False
811 self.byProfiles = False
812 # self.byProfiles = False
812
813 #
813 self.bloques = None
814 # self.bloques = None
814 self.bloque0 = None
815 # self.bloque0 = None
815
816 #
816 self.index = 0
817 # self.index = 0
817
818 #
818 self.buffer = 0
819 # self.buffer = 0
819 self.buffer2 = 0
820 # self.buffer2 = 0
820 self.buffer3 = 0
821 # self.buffer3 = 0
821
822 #
822
823 #
823 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
824 # def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
824
825 #
825 self.nChannels = dataOut.nChannels
826 # self.nChannels = dataOut.nChannels
826 self.nProf = dataOut.nProfiles
827 # self.nProf = dataOut.nProfiles
827 self.nPairs = dataOut.data_cspc.shape[0]
828 # self.nPairs = dataOut.data_cspc.shape[0]
828 self.pairsArray = numpy.array(dataOut.pairsList)
829 # self.pairsArray = numpy.array(dataOut.pairsList)
829 self.spectra = dataOut.data_spc
830 # self.spectra = dataOut.data_spc
830 self.cspectra = dataOut.data_cspc
831 # self.cspectra = dataOut.data_cspc
831 self.heights = dataOut.heightList #alturas totales
832 # self.heights = dataOut.heightList #alturas totales
832 self.nHeights = len(self.heights)
833 # self.nHeights = len(self.heights)
833 self.min_hei = min_hei
834 # self.min_hei = min_hei
834 self.max_hei = max_hei
835 # self.max_hei = max_hei
835 if (self.min_hei == None):
836 # if (self.min_hei == None):
836 self.min_hei = 0
837 # self.min_hei = 0
837 if (self.max_hei == None):
838 # if (self.max_hei == None):
838 self.max_hei = dataOut.heightList[-1]
839 # self.max_hei = dataOut.heightList[-1]
839 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
840 # self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
840 self.heightsClean = self.heights[self.hval] #alturas filtradas
841 # self.heightsClean = self.heights[self.hval] #alturas filtradas
841 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
842 # self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
842 self.nHeightsClean = len(self.heightsClean)
843 # self.nHeightsClean = len(self.heightsClean)
843 self.channels = dataOut.channelList
844 # self.channels = dataOut.channelList
844 self.nChan = len(self.channels)
845 # self.nChan = len(self.channels)
845 self.nIncohInt = dataOut.nIncohInt
846 # self.nIncohInt = dataOut.nIncohInt
846 self.__initime = dataOut.utctime
847 # self.__initime = dataOut.utctime
847 self.maxAltInd = self.hval[-1]+1
848 # self.maxAltInd = self.hval[-1]+1
848 self.minAltInd = self.hval[0]
849 # self.minAltInd = self.hval[0]
849
850 #
850 self.crosspairs = dataOut.pairsList
851 # self.crosspairs = dataOut.pairsList
851 self.nPairs = len(self.crosspairs)
852 # self.nPairs = len(self.crosspairs)
852 self.normFactor = dataOut.normFactor
853 # self.normFactor = dataOut.normFactor
853 self.nFFTPoints = dataOut.nFFTPoints
854 # self.nFFTPoints = dataOut.nFFTPoints
854 self.ippSeconds = dataOut.ippSeconds
855 # self.ippSeconds = dataOut.ippSeconds
855 self.currentTime = self.__initime
856 # self.currentTime = self.__initime
856 self.pairsArray = numpy.array(dataOut.pairsList)
857 # self.pairsArray = numpy.array(dataOut.pairsList)
857 self.factor_stdv = factor_stdv
858 # self.factor_stdv = factor_stdv
858
859 #
859 if n != None :
860 # if n != None :
860 self.byProfiles = True
861 # self.byProfiles = True
861 self.nIntProfiles = n
862 # self.nIntProfiles = n
862 else:
863 # else:
863 self.__integrationtime = timeInterval
864 # self.__integrationtime = timeInterval
864
865 #
865 self.__dataReady = False
866 # self.__dataReady = False
866 self.isConfig = True
867 # self.isConfig = True
867
868 #
868
869 #
869
870 #
870 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
871 # def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
871 #print("runing cleanRayleigh")
872 # #print("runing cleanRayleigh")
872 if not self.isConfig :
873 # if not self.isConfig :
873
874 #
874 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
875 # self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
875
876 #
876 tini=dataOut.utctime
877 # tini=dataOut.utctime
877
878 #
878 if self.byProfiles:
879 # if self.byProfiles:
879 if self.__profIndex == self.nIntProfiles:
880 # if self.__profIndex == self.nIntProfiles:
880 self.__dataReady = True
881 # self.__dataReady = True
881 else:
882 # else:
882 if (tini - self.__initime) >= self.__integrationtime:
883 # if (tini - self.__initime) >= self.__integrationtime:
883
884 #
884 self.__dataReady = True
885 # self.__dataReady = True
885 self.__initime = tini
886 # self.__initime = tini
886
887 #
887 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
888 # #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
888
889 #
889 if self.__dataReady:
890 # if self.__dataReady:
890
891 #
891 self.__profIndex = 0
892 # self.__profIndex = 0
892 jspc = self.buffer
893 # jspc = self.buffer
893 jcspc = self.buffer2
894 # jcspc = self.buffer2
894 #jnoise = self.buffer3
895 # #jnoise = self.buffer3
895 self.buffer = dataOut.data_spc
896 # self.buffer = dataOut.data_spc
896 self.buffer2 = dataOut.data_cspc
897 # self.buffer2 = dataOut.data_cspc
897 #self.buffer3 = dataOut.noise
898 # #self.buffer3 = dataOut.noise
898 self.currentTime = dataOut.utctime
899 # self.currentTime = dataOut.utctime
899 if numpy.any(jspc) :
900 # if numpy.any(jspc) :
900 #print( jspc.shape, jcspc.shape)
901 # #print( jspc.shape, jcspc.shape)
901 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
902 # jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
902 try:
903 # try:
903 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
904 # jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
904 except:
905 # except:
905 print("no cspc")
906 # print("no cspc")
906 self.__dataReady = False
907 # self.__dataReady = False
907 #print( jspc.shape, jcspc.shape)
908 # #print( jspc.shape, jcspc.shape)
908 dataOut.flagNoData = False
909 # dataOut.flagNoData = False
909 else:
910 # else:
910 dataOut.flagNoData = True
911 # dataOut.flagNoData = True
911 self.__dataReady = False
912 # self.__dataReady = False
912 return dataOut
913 # return dataOut
913 else:
914 # else:
914 #print( len(self.buffer))
915 # #print( len(self.buffer))
915 if numpy.any(self.buffer):
916 # if numpy.any(self.buffer):
916 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
917 # self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
917 try:
918 # try:
918 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
919 # self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
919 self.buffer3 += dataOut.data_dc
920 # self.buffer3 += dataOut.data_dc
920 except:
921 # except:
921 pass
922 # pass
922 else:
923 # else:
923 self.buffer = dataOut.data_spc
924 # self.buffer = dataOut.data_spc
924 self.buffer2 = dataOut.data_cspc
925 # self.buffer2 = dataOut.data_cspc
925 self.buffer3 = dataOut.data_dc
926 # self.buffer3 = dataOut.data_dc
926 #print self.index, self.fint
927 # #print self.index, self.fint
927 #print self.buffer2.shape
928 # #print self.buffer2.shape
928 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
929 # dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
929 self.__profIndex += 1
930 # self.__profIndex += 1
930 return dataOut ## NOTE: REV
931 # return dataOut ## NOTE: REV
931
932 #
932
933 #
933 #index = tini.tm_hour*12+tini.tm_min/5
934 # #index = tini.tm_hour*12+tini.tm_min/5
934 '''
935 # '''
935 #REVISAR
936 # #REVISAR
936 '''
937 # '''
937 # jspc = jspc/self.nFFTPoints/self.normFactor
938 # # jspc = jspc/self.nFFTPoints/self.normFactor
938 # jcspc = jcspc/self.nFFTPoints/self.normFactor
939 # # jcspc = jcspc/self.nFFTPoints/self.normFactor
939
940 #
940
941 #
941
942 #
942 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
943 # tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
943 dataOut.data_spc = tmp_spectra
944 # dataOut.data_spc = tmp_spectra
944 dataOut.data_cspc = tmp_cspectra
945 # dataOut.data_cspc = tmp_cspectra
945
946 #
946 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
947 # #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
947
948 #
948 dataOut.data_dc = self.buffer3
949 # dataOut.data_dc = self.buffer3
949 dataOut.nIncohInt *= self.nIntProfiles
950 # dataOut.nIncohInt *= self.nIntProfiles
950 dataOut.max_nIncohInt = self.nIntProfiles
951 # dataOut.max_nIncohInt = self.nIntProfiles
951 dataOut.utctime = self.currentTime #tiempo promediado
952 # dataOut.utctime = self.currentTime #tiempo promediado
952 #print("Time: ",time.localtime(dataOut.utctime))
953 # #print("Time: ",time.localtime(dataOut.utctime))
953 # dataOut.data_spc = sat_spectra
954 # # dataOut.data_spc = sat_spectra
954 # dataOut.data_cspc = sat_cspectra
955 # # dataOut.data_cspc = sat_cspectra
955 self.buffer = 0
956 # self.buffer = 0
956 self.buffer2 = 0
957 # self.buffer2 = 0
957 self.buffer3 = 0
958 # self.buffer3 = 0
958
959 #
959 return dataOut
960 # return dataOut
960
961 #
961 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
962 # def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
962 print("OP cleanRayleigh")
963 # print("OP cleanRayleigh")
963 #import matplotlib.pyplot as plt
964 # #import matplotlib.pyplot as plt
964 #for k in range(149):
965 # #for k in range(149):
965 #channelsProcssd = []
966 # #channelsProcssd = []
966 #channelA_ok = False
967 # #channelA_ok = False
967 #rfunc = cspectra.copy() #self.bloques
968 # #rfunc = cspectra.copy() #self.bloques
968 rfunc = spectra.copy()
969 # rfunc = spectra.copy()
969 #rfunc = cspectra
970 # #rfunc = cspectra
970 #val_spc = spectra*0.0 #self.bloque0*0.0
971 # #val_spc = spectra*0.0 #self.bloque0*0.0
971 #val_cspc = cspectra*0.0 #self.bloques*0.0
972 # #val_cspc = cspectra*0.0 #self.bloques*0.0
972 #in_sat_spectra = spectra.copy() #self.bloque0
973 # #in_sat_spectra = spectra.copy() #self.bloque0
973 #in_sat_cspectra = cspectra.copy() #self.bloques
974 # #in_sat_cspectra = cspectra.copy() #self.bloques
974
975 #
975
976 #
976 ###ONLY FOR TEST:
977 # ###ONLY FOR TEST:
977 raxs = math.ceil(math.sqrt(self.nPairs))
978 # raxs = math.ceil(math.sqrt(self.nPairs))
978 if raxs == 0:
979 # if raxs == 0:
979 raxs = 1
980 # raxs = 1
980 caxs = math.ceil(self.nPairs/raxs)
981 # caxs = math.ceil(self.nPairs/raxs)
981 if self.nPairs <4:
982 # if self.nPairs <4:
982 raxs = 2
983 # raxs = 2
983 caxs = 2
984 # caxs = 2
984 #print(raxs, caxs)
985 # #print(raxs, caxs)
985 fft_rev = 14 #nFFT to plot
986 # fft_rev = 14 #nFFT to plot
986 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
987 # hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
987 hei_rev = hei_rev[0]
988 # hei_rev = hei_rev[0]
988 #print(hei_rev)
989 # #print(hei_rev)
989
990 #
990 #print numpy.absolute(rfunc[:,0,0,14])
991 # #print numpy.absolute(rfunc[:,0,0,14])
991
992 #
992 gauss_fit, covariance = None, None
993 # gauss_fit, covariance = None, None
993 for ih in range(self.minAltInd,self.maxAltInd):
994 # for ih in range(self.minAltInd,self.maxAltInd):
994 for ifreq in range(self.nFFTPoints):
995 # for ifreq in range(self.nFFTPoints):
995 '''
996 # '''
996 ###ONLY FOR TEST:
997 # ###ONLY FOR TEST:
997 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
998 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
998 fig, axs = plt.subplots(raxs, caxs)
999 # fig, axs = plt.subplots(raxs, caxs)
999 fig2, axs2 = plt.subplots(raxs, caxs)
1000 # fig2, axs2 = plt.subplots(raxs, caxs)
1000 col_ax = 0
1001 # col_ax = 0
1001 row_ax = 0
1002 # row_ax = 0
1002 '''
1003 # '''
1003 #print(self.nPairs)
1004 # #print(self.nPairs)
1004 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
1005 # for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
1005 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
1006 # # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
1006 # continue
1007 # # continue
1007 # if not self.crosspairs[ii][0] in channelsProcssd:
1008 # # if not self.crosspairs[ii][0] in channelsProcssd:
1008 # channelA_ok = True
1009 # # channelA_ok = True
1009 #print("pair: ",self.crosspairs[ii])
1010 # #print("pair: ",self.crosspairs[ii])
1010 '''
1011 # '''
1011 ###ONLY FOR TEST:
1012 # ###ONLY FOR TEST:
1012 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
1013 # if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
1013 col_ax = 0
1014 # col_ax = 0
1014 row_ax += 1
1015 # row_ax += 1
1015 '''
1016 # '''
1016 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
1017 # func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
1017 #print(func2clean.shape)
1018 # #print(func2clean.shape)
1018 val = (numpy.isfinite(func2clean)==True).nonzero()
1019 # val = (numpy.isfinite(func2clean)==True).nonzero()
1019
1020 #
1020 if len(val)>0: #limitador
1021 # if len(val)>0: #limitador
1021 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
1022 # min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
1022 if min_val <= -40 :
1023 # if min_val <= -40 :
1023 min_val = -40
1024 # min_val = -40
1024 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
1025 # max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
1025 if max_val >= 200 :
1026 # if max_val >= 200 :
1026 max_val = 200
1027 # max_val = 200
1027 #print min_val, max_val
1028 # #print min_val, max_val
1028 step = 1
1029 # step = 1
1029 #print("Getting bins and the histogram")
1030 # #print("Getting bins and the histogram")
1030 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
1031 # x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
1031 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1032 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1032 #print(len(y_dist),len(binstep[:-1]))
1033 # #print(len(y_dist),len(binstep[:-1]))
1033 #print(row_ax,col_ax, " ..")
1034 # #print(row_ax,col_ax, " ..")
1034 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
1035 # #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
1035 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
1036 # mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
1036 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
1037 # sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
1037 parg = [numpy.amax(y_dist),mean,sigma]
1038 # parg = [numpy.amax(y_dist),mean,sigma]
1038
1039 #
1039 newY = None
1040 # newY = None
1040
1041 #
1041 try :
1042 # try :
1042 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
1043 # gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
1043 mode = gauss_fit[1]
1044 # mode = gauss_fit[1]
1044 stdv = gauss_fit[2]
1045 # stdv = gauss_fit[2]
1045 #print(" FIT OK",gauss_fit)
1046 # #print(" FIT OK",gauss_fit)
1046 '''
1047 # '''
1047 ###ONLY FOR TEST:
1048 # ###ONLY FOR TEST:
1048 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1049 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1049 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
1050 # newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
1050 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1051 # axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1051 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1052 # axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1052 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1053 # axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1053 '''
1054 # '''
1054 except:
1055 # except:
1055 mode = mean
1056 # mode = mean
1056 stdv = sigma
1057 # stdv = sigma
1057 #print("FIT FAIL")
1058 # #print("FIT FAIL")
1058 #continue
1059 # #continue
1059
1060 #
1060
1061 #
1061 #print(mode,stdv)
1062 # #print(mode,stdv)
1062 #Removing echoes greater than mode + std_factor*stdv
1063 # #Removing echoes greater than mode + std_factor*stdv
1063 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
1064 # noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
1064 #noval tiene los indices que se van a remover
1065 # #noval tiene los indices que se van a remover
1065 #print("Chan ",ii," novals: ",len(noval[0]))
1066 # #print("Chan ",ii," novals: ",len(noval[0]))
1066 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
1067 # if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
1067 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
1068 # novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
1068 #print(novall)
1069 # #print(novall)
1069 #print(" ",self.pairsArray[ii])
1070 # #print(" ",self.pairsArray[ii])
1070 #cross_pairs = self.pairsArray[ii]
1071 # #cross_pairs = self.pairsArray[ii]
1071 #Getting coherent echoes which are removed.
1072 # #Getting coherent echoes which are removed.
1072 # if len(novall[0]) > 0:
1073 # # if len(novall[0]) > 0:
1073 #
1074 # #
1074 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
1075 # # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
1075 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
1076 # # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
1076 # val_cspc[novall[0],ii,ifreq,ih] = 1
1077 # # val_cspc[novall[0],ii,ifreq,ih] = 1
1077 #print("OUT NOVALL 1")
1078 # #print("OUT NOVALL 1")
1078 try:
1079 # try:
1079 pair = (self.channels[ii],self.channels[ii + 1])
1080 # pair = (self.channels[ii],self.channels[ii + 1])
1080 except:
1081 # except:
1081 pair = (99,99)
1082 # pair = (99,99)
1082 #print("par ", pair)
1083 # #print("par ", pair)
1083 if ( pair in self.crosspairs):
1084 # if ( pair in self.crosspairs):
1084 q = self.crosspairs.index(pair)
1085 # q = self.crosspairs.index(pair)
1085 #print("estΓ‘ aqui: ", q, (ii,ii + 1))
1086 # #print("estΓ‘ aqui: ", q, (ii,ii + 1))
1086 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
1087 # new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
1087 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
1088 # cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
1088
1089 #
1089 #if channelA_ok:
1090 # #if channelA_ok:
1090 #chA = self.channels.index(cross_pairs[0])
1091 # #chA = self.channels.index(cross_pairs[0])
1091 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
1092 # new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
1092 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
1093 # spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
1093 #channelA_ok = False
1094 # #channelA_ok = False
1094
1095 #
1095 # chB = self.channels.index(cross_pairs[1])
1096 # # chB = self.channels.index(cross_pairs[1])
1096 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
1097 # # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
1097 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
1098 # # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
1098 #
1099 # #
1099 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
1100 # # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
1100 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
1101 # # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
1101 '''
1102 # '''
1102 ###ONLY FOR TEST:
1103 # ###ONLY FOR TEST:
1103 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1104 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1104 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
1105 # func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
1105 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1106 # y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
1106 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1107 # axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
1107 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1108 # axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
1108 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1109 # axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
1109 '''
1110 # '''
1110 '''
1111 # '''
1111 ###ONLY FOR TEST:
1112 # ###ONLY FOR TEST:
1112 col_ax += 1 #contador de ploteo columnas
1113 # col_ax += 1 #contador de ploteo columnas
1113 ##print(col_ax)
1114 # ##print(col_ax)
1114 ###ONLY FOR TEST:
1115 # ###ONLY FOR TEST:
1115 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1116 # if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
1116 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
1117 # title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
1117 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
1118 # title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
1118 fig.suptitle(title)
1119 # fig.suptitle(title)
1119 fig2.suptitle(title2)
1120 # fig2.suptitle(title2)
1120 plt.show()
1121 # plt.show()
1121 '''
1122 # '''
1122 ##################################################################################################
1123 # ##################################################################################################
1123
1124 #
1124 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
1125 # #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
1125 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
1126 # out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
1126 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
1127 # out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
1127 for ih in range(self.nHeights):
1128 # for ih in range(self.nHeights):
1128 for ifreq in range(self.nFFTPoints):
1129 # for ifreq in range(self.nFFTPoints):
1129 for ich in range(self.nChan):
1130 # for ich in range(self.nChan):
1130 tmp = spectra[:,ich,ifreq,ih]
1131 # tmp = spectra[:,ich,ifreq,ih]
1131 valid = (numpy.isfinite(tmp[:])==True).nonzero()
1132 # valid = (numpy.isfinite(tmp[:])==True).nonzero()
1132
1133 #
1133 if len(valid[0]) >0 :
1134 # if len(valid[0]) >0 :
1134 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1135 # out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1135
1136 #
1136 for icr in range(self.nPairs):
1137 # for icr in range(self.nPairs):
1137 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
1138 # tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
1138 valid = (numpy.isfinite(tmp)==True).nonzero()
1139 # valid = (numpy.isfinite(tmp)==True).nonzero()
1139 if len(valid[0]) > 0:
1140 # if len(valid[0]) > 0:
1140 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1141 # out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
1141
1142 #
1142 return out_spectra, out_cspectra
1143 # return out_spectra, out_cspectra
1143
1144 #
1144 def REM_ISOLATED_POINTS(self,array,rth):
1145 # def REM_ISOLATED_POINTS(self,array,rth):
1145 # import matplotlib.pyplot as plt
1146 # # import matplotlib.pyplot as plt
1146 if rth == None :
1147 # if rth == None :
1147 rth = 4
1148 # rth = 4
1148 #print("REM ISO")
1149 # #print("REM ISO")
1149 num_prof = len(array[0,:,0])
1150 # num_prof = len(array[0,:,0])
1150 num_hei = len(array[0,0,:])
1151 # num_hei = len(array[0,0,:])
1151 n2d = len(array[:,0,0])
1152 # n2d = len(array[:,0,0])
1152
1153 #
1153 for ii in range(n2d) :
1154 # for ii in range(n2d) :
1154 #print ii,n2d
1155 # #print ii,n2d
1155 tmp = array[ii,:,:]
1156 # tmp = array[ii,:,:]
1156 #print tmp.shape, array[ii,101,:],array[ii,102,:]
1157 # #print tmp.shape, array[ii,101,:],array[ii,102,:]
1157
1158 #
1158 # fig = plt.figure(figsize=(6,5))
1159 # # fig = plt.figure(figsize=(6,5))
1159 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1160 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1160 # ax = fig.add_axes([left, bottom, width, height])
1161 # # ax = fig.add_axes([left, bottom, width, height])
1161 # x = range(num_prof)
1162 # # x = range(num_prof)
1162 # y = range(num_hei)
1163 # # y = range(num_hei)
1163 # cp = ax.contour(y,x,tmp)
1164 # # cp = ax.contour(y,x,tmp)
1164 # ax.clabel(cp, inline=True,fontsize=10)
1165 # # ax.clabel(cp, inline=True,fontsize=10)
1165 # plt.show()
1166 # # plt.show()
1166
1167 #
1167 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
1168 # #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
1168 tmp = numpy.reshape(tmp,num_prof*num_hei)
1169 # tmp = numpy.reshape(tmp,num_prof*num_hei)
1169 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
1170 # indxs1 = (numpy.isfinite(tmp)==True).nonzero()
1170 indxs2 = (tmp > 0).nonzero()
1171 # indxs2 = (tmp > 0).nonzero()
1171
1172 #
1172 indxs1 = (indxs1[0])
1173 # indxs1 = (indxs1[0])
1173 indxs2 = indxs2[0]
1174 # indxs2 = indxs2[0]
1174 #indxs1 = numpy.array(indxs1[0])
1175 # #indxs1 = numpy.array(indxs1[0])
1175 #indxs2 = numpy.array(indxs2[0])
1176 # #indxs2 = numpy.array(indxs2[0])
1176 indxs = None
1177 # indxs = None
1177 #print indxs1 , indxs2
1178 # #print indxs1 , indxs2
1178 for iv in range(len(indxs2)):
1179 # for iv in range(len(indxs2)):
1179 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
1180 # indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
1180 #print len(indxs2), indv
1181 # #print len(indxs2), indv
1181 if len(indv[0]) > 0 :
1182 # if len(indv[0]) > 0 :
1182 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
1183 # indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
1183 # print indxs
1184 # # print indxs
1184 indxs = indxs[1:]
1185 # indxs = indxs[1:]
1185 #print(indxs, len(indxs))
1186 # #print(indxs, len(indxs))
1186 if len(indxs) < 4 :
1187 # if len(indxs) < 4 :
1187 array[ii,:,:] = 0.
1188 # array[ii,:,:] = 0.
1188 return
1189 # return
1189
1190 #
1190 xpos = numpy.mod(indxs ,num_hei)
1191 # xpos = numpy.mod(indxs ,num_hei)
1191 ypos = (indxs / num_hei)
1192 # ypos = (indxs / num_hei)
1192 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1193 # sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
1193 #print sx
1194 # #print sx
1194 xpos = xpos[sx]
1195 # xpos = xpos[sx]
1195 ypos = ypos[sx]
1196 # ypos = ypos[sx]
1196
1197 #
1197 # *********************************** Cleaning isolated points **********************************
1198 # # *********************************** Cleaning isolated points **********************************
1198 ic = 0
1199 # ic = 0
1199 while True :
1200 # while True :
1200 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1201 # r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
1201 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1202 # #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
1202 #plt.plot(r)
1203 # #plt.plot(r)
1203 #plt.show()
1204 # #plt.show()
1204 no_coh1 = (numpy.isfinite(r)==True).nonzero()
1205 # no_coh1 = (numpy.isfinite(r)==True).nonzero()
1205 no_coh2 = (r <= rth).nonzero()
1206 # no_coh2 = (r <= rth).nonzero()
1206 #print r, no_coh1, no_coh2
1207 # #print r, no_coh1, no_coh2
1207 no_coh1 = numpy.array(no_coh1[0])
1208 # no_coh1 = numpy.array(no_coh1[0])
1208 no_coh2 = numpy.array(no_coh2[0])
1209 # no_coh2 = numpy.array(no_coh2[0])
1209 no_coh = None
1210 # no_coh = None
1210 #print valid1 , valid2
1211 # #print valid1 , valid2
1211 for iv in range(len(no_coh2)):
1212 # for iv in range(len(no_coh2)):
1212 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1213 # indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
1213 if len(indv[0]) > 0 :
1214 # if len(indv[0]) > 0 :
1214 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1215 # no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
1215 no_coh = no_coh[1:]
1216 # no_coh = no_coh[1:]
1216 #print len(no_coh), no_coh
1217 # #print len(no_coh), no_coh
1217 if len(no_coh) < 4 :
1218 # if len(no_coh) < 4 :
1218 #print xpos[ic], ypos[ic], ic
1219 # #print xpos[ic], ypos[ic], ic
1219 # plt.plot(r)
1220 # # plt.plot(r)
1220 # plt.show()
1221 # # plt.show()
1221 xpos[ic] = numpy.nan
1222 # xpos[ic] = numpy.nan
1222 ypos[ic] = numpy.nan
1223 # ypos[ic] = numpy.nan
1223
1224 #
1224 ic = ic + 1
1225 # ic = ic + 1
1225 if (ic == len(indxs)) :
1226 # if (ic == len(indxs)) :
1226 break
1227 # break
1227 #print( xpos, ypos)
1228 # #print( xpos, ypos)
1228
1229 #
1229 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1230 # indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1230 #print indxs[0]
1231 # #print indxs[0]
1231 if len(indxs[0]) < 4 :
1232 # if len(indxs[0]) < 4 :
1232 array[ii,:,:] = 0.
1233 # array[ii,:,:] = 0.
1233 return
1234 # return
1234
1235 #
1235 xpos = xpos[indxs[0]]
1236 # xpos = xpos[indxs[0]]
1236 ypos = ypos[indxs[0]]
1237 # ypos = ypos[indxs[0]]
1237 for i in range(0,len(ypos)):
1238 # for i in range(0,len(ypos)):
1238 ypos[i]=int(ypos[i])
1239 # ypos[i]=int(ypos[i])
1239 junk = tmp
1240 # junk = tmp
1240 tmp = junk*0.0
1241 # tmp = junk*0.0
1241
1242 #
1242 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1243 # tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1243 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1244 # array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1244
1245 #
1245 #print array.shape
1246 # #print array.shape
1246 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1247 # #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1247 #print tmp.shape
1248 # #print tmp.shape
1248
1249 #
1249 # fig = plt.figure(figsize=(6,5))
1250 # # fig = plt.figure(figsize=(6,5))
1250 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1251 # # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1251 # ax = fig.add_axes([left, bottom, width, height])
1252 # # ax = fig.add_axes([left, bottom, width, height])
1252 # x = range(num_prof)
1253 # # x = range(num_prof)
1253 # y = range(num_hei)
1254 # # y = range(num_hei)
1254 # cp = ax.contour(y,x,array[ii,:,:])
1255 # # cp = ax.contour(y,x,array[ii,:,:])
1255 # ax.clabel(cp, inline=True,fontsize=10)
1256 # # ax.clabel(cp, inline=True,fontsize=10)
1256 # plt.show()
1257 # # plt.show()
1257 return array
1258 # return array
1258
1259 #
1259
1260
1260 class IntegrationFaradaySpectra(Operation):
1261 class IntegrationFaradaySpectra(Operation):
1261
1262
@@ -1313,7 +1314,7 class IntegrationFaradaySpectra(Operation):
1313 self.navg = avg
1314 self.navg = avg
1314 #self.ByLags = dataOut.ByLags ###REDEFINIR
1315 #self.ByLags = dataOut.ByLags ###REDEFINIR
1315 self.ByLags = False
1316 self.ByLags = False
1316 self.maxProfilesInt = 1
1317 self.maxProfilesInt = 0
1317 self.__nChannels = dataOut.nChannels
1318 self.__nChannels = dataOut.nChannels
1318 if DPL != None:
1319 if DPL != None:
1319 self.DPL=DPL
1320 self.DPL=DPL
@@ -1542,7 +1543,7 class IntegrationFaradaySpectra(Operation):
1542 data_cspc = None
1543 data_cspc = None
1543 data_dc = self.__buffer_dc
1544 data_dc = self.__buffer_dc
1544 #(CH, HEIGH)
1545 #(CH, HEIGH)
1545 self.maxProfilesInt = self.__profIndex
1546 self.maxProfilesInt = self.__profIndex - 1
1546 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1547 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1547
1548
1548 self.__buffer_spc = []
1549 self.__buffer_spc = []
@@ -1656,8 +1657,8 class IntegrationFaradaySpectra(Operation):
1656
1657
1657
1658
1658 self.dataOut.nIncohInt *= self.n_ints
1659 self.dataOut.nIncohInt *= self.n_ints
1659 self.dataOut.max_nIncohInt = self.maxProfilesInt
1660 #print("maxProfilesInt: ",self.maxProfilesInt)
1660 #print(self.dataOut.max_nIncohInt)
1661
1661 self.dataOut.utctime = avgdatatime
1662 self.dataOut.utctime = avgdatatime
1662 self.dataOut.flagNoData = False
1663 self.dataOut.flagNoData = False
1663 #print("Faraday Integration DONE...", self.dataOut.data_cspc)
1664 #print("Faraday Integration DONE...", self.dataOut.data_cspc)
@@ -2130,7 +2131,6 class IncohInt(Operation):
2130 dataOut.data_outlier = self.nOutliers
2131 dataOut.data_outlier = self.nOutliers
2131 dataOut.utctime = avgdatatime
2132 dataOut.utctime = avgdatatime
2132 dataOut.flagNoData = False
2133 dataOut.flagNoData = False
2133 dataOut.max_nIncohInt += self.__profIndex
2134 self.incohInt = 0
2134 self.incohInt = 0
2135 self.nOutliers = 0
2135 self.nOutliers = 0
2136 self.__profIndex = 0
2136 self.__profIndex = 0
Show file before | Show file after | Hide comments
@@ -3250,7 +3250,7 class RemoveProfileSats2(Operation):
3250 fnoise = scipy.signal.filtfilt(b, a, noise_ref)
3250 fnoise = scipy.signal.filtfilt(b, a, noise_ref)
3251 #noise_refdB = 10* numpy.log10(noise_ref)
3251 #noise_refdB = 10* numpy.log10(noise_ref)
3252 #print("Noise ",numpy.percentile(noise_ref,95))
3252 #print("Noise ",numpy.percentile(noise_ref,95))
3253 p85 = numpy.percentile(fnoise,85)
3253 p85 = numpy.percentile(fnoise,83)
3254 mean_noise = fnoise.mean()
3254 mean_noise = fnoise.mean()
3255 if self.prev_pnoise != None:
3255 if self.prev_pnoise != None:
3256 if mean_noise < (1.5 * self.prev_pnoise) :
3256 if mean_noise < (1.5 * self.prev_pnoise) :
@@ -3280,7 +3280,7 class RemoveProfileSats2(Operation):
3280 if index[0].size < int(self.navg*profiles):
3280 if index[0].size < int(self.navg*profiles):
3281 indexes = numpy.append(indexes, index[0])
3281 indexes = numpy.append(indexes, index[0])
3282
3282
3283 #print(indexes)
3283
3284 # from matplotlib import pyplot as plt
3284 # from matplotlib import pyplot as plt
3285 # fig, ax = plt.subplots()
3285 # fig, ax = plt.subplots()
3286 # ax.plot(fpower)
3286 # ax.plot(fpower)
@@ -3288,7 +3288,8 class RemoveProfileSats2(Operation):
3288 # ax.axhline(std, color='b')
3288 # ax.axhline(std, color='b')
3289 # plt.grid()
3289 # plt.grid()
3290 # plt.show()
3290 # plt.show()
3291 # print(indexes)
3291
3292 #print(indexes)
3292
3293
3293 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3294 #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
3294 #outliers_IDs = numpy.unique(outliers_IDs)
3295 #outliers_IDs = numpy.unique(outliers_IDs)
@@ -3301,15 +3302,19 class RemoveProfileSats2(Operation):
3301
3302
3302
3303
3303 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
3304 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
3304 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
3305 hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier) #es outlier
3305 hist_outliers_indexes = hist_outliers_indexes[0]
3306 hist_outliers_indexes = hist_outliers_indexes[0]
3306 if len(hist_outliers_indexes>0):
3307 if len(hist_outliers_indexes>0):
3307 hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
3308 hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
3308
3309
3309 bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
3310 bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
3310 outlier_loc_index = []
3311 outlier_loc_index = []
3311
3312 #print("out indexes ", bins_outliers_indexes)
3312 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] + self.profileMargin) ]
3313 if len(bins_outliers_indexes) <= 3:
3314 extprof = 0
3315 else:
3316 extprof = self.profileMargin
3317 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) ]
3313 outlier_loc_index = numpy.asarray(outlier_loc_index)
3318 outlier_loc_index = numpy.asarray(outlier_loc_index)
3314 # if len(outlier_loc_index)>1:
3319 # if len(outlier_loc_index)>1:
3315 # ipmax = numpy.where(fpower==fpower.max())[0]
3320 # ipmax = numpy.where(fpower==fpower.max())[0]
@@ -3331,7 +3336,7 class RemoveProfileSats2(Operation):
3331 # o = numpy.nanstd(dat)
3336 # o = numpy.nanstd(dat)
3332 # #print(m, o, x.shape, y.shape)
3337 # #print(m, o, x.shape, y.shape)
3333 # #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
3338 # #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
3334 # c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = 55, vmax = 65)
3339 # c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70)
3335 # ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3340 # ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w')
3336 # #fig.colorbar(c)
3341 # #fig.colorbar(c)
3337 # ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
3342 # ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')
General Comments 0
You need to be logged in to leave comments. Login now