| @@ -9,7 +9,7 | |||
|
|
9 | 9 | import numpy |
|
|
10 | 10 | |
|
|
11 | 11 | from schainpy.model.graphics.jroplot_base import Plot, plt |
|
|
12 | ||
|
|
12 | import matplotlib.pyplot as plt | |
|
|
13 | 13 | |
|
|
14 | 14 | class SpectraHeisPlot(Plot): |
|
|
15 | 15 | |
| @@ -39,11 +39,12 class SpectraHeisPlot(Plot): | |||
|
|
39 | 39 | def plot(self): |
|
|
40 | 40 | |
|
|
41 | 41 | c = 3E8 |
|
|
42 | deltaHeight = self.data.yrange[1] - self.data.yrange[0] | |
|
|
42 | deltaHeight = self.data.yrange[1] - self.data.yrange[0] # yrange = heightList | |
|
|
43 | 43 | x = numpy.arange(-1*len(self.data.yrange)/2., len(self.data.yrange)/2.)*(c/(2*deltaHeight*len(self.data.yrange)*1000)) |
|
|
44 | 44 | self.y = self.data[-1]['spc_heis'] |
|
|
45 | 45 | self.titles = [] |
|
|
46 | 46 | |
|
|
47 | Maintitle = "Range from %d km to %d km" %(int(self.data.yrange[0]),int(self.data.yrange[-1])) | |
|
|
47 | 48 | for n, ax in enumerate(self.axes): |
|
|
48 | 49 | ychannel = self.y[n,:] |
|
|
49 | 50 | if ax.firsttime: |
| @@ -54,7 +55,7 class SpectraHeisPlot(Plot): | |||
|
|
54 | 55 | ax.plt.set_data(x, ychannel) |
|
|
55 | 56 | |
|
|
56 | 57 | self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel))) |
|
|
57 | ||
|
|
58 | plt.suptitle(Maintitle) | |
|
|
58 | 59 | |
|
|
59 | 60 | class RTIHeisPlot(Plot): |
|
|
60 | 61 | |
| @@ -22,11 +22,11 class ScopePlot(Plot): | |||
|
|
22 | 22 | def setup(self): |
|
|
23 | 23 | |
|
|
24 | 24 | self.xaxis = 'Range (Km)' |
|
|
25 | self.ncols = 1 | |
|
|
26 | self.nrows = 1 | |
|
|
27 | self.nplots = 1 | |
|
|
25 | self.nplots = len(self.data.channels) | |
|
|
26 | self.nrows = int(numpy.ceil(self.nplots/2)) | |
|
|
27 | self.ncols = int(numpy.ceil(self.nplots/self.nrows)) | |
|
|
28 | 28 | self.ylabel = 'Intensity [dB]' |
|
|
29 | self.titles = ['Scope'] | |
|
|
29 | self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels] | |
|
|
30 | 30 | self.colorbar = False |
|
|
31 | 31 | self.width = 6 |
|
|
32 | 32 | self.height = 4 |
| @@ -56,15 +56,13 class ScopePlot(Plot): | |||
|
|
56 | 56 | |
|
|
57 | 57 | yreal = y[channelIndexList,:].real |
|
|
58 | 58 | yimag = y[channelIndexList,:].imag |
|
|
59 | title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) | |
|
|
59 | Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) | |
|
|
60 | 60 | self.xlabel = "Range (Km)" |
|
|
61 | 61 | self.ylabel = "Intensity - IQ" |
|
|
62 | 62 | |
|
|
63 | 63 | self.y = yreal |
|
|
64 | 64 | self.x = x |
|
|
65 | 65 | |
|
|
66 | self.titles[0] = title | |
|
|
67 | ||
|
|
68 | 66 | for i,ax in enumerate(self.axes): |
|
|
69 | 67 | title = "Channel %d" %(i) |
|
|
70 | 68 | if ax.firsttime: |
| @@ -75,6 +73,7 class ScopePlot(Plot): | |||
|
|
75 | 73 | else: |
|
|
76 | 74 | ax.plt_r.set_data(x, yreal[i,:]) |
|
|
77 | 75 | ax.plt_i.set_data(x, yimag[i,:]) |
|
|
76 | plt.suptitle(Maintitle) | |
|
|
78 | 77 | |
|
|
79 | 78 | def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle): |
|
|
80 | 79 | y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) |
| @@ -99,6 +98,7 class ScopePlot(Plot): | |||
|
|
99 | 98 | else: |
|
|
100 | 99 | ax.plt_r.set_data(x, ychannel) |
|
|
101 | 100 | |
|
|
101 | ||
|
|
102 | 102 | def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle): |
|
|
103 | 103 | |
|
|
104 | 104 | |
| @@ -689,9 +689,9 class Reader(object): | |||
|
|
689 | 689 | @staticmethod |
|
|
690 | 690 | def isDateTimeInRange(dt, startDate, endDate, startTime, endTime): |
|
|
691 | 691 | """Check if the given datetime is in range""" |
|
|
692 | ||
|
|
693 | if startDate <= dt.date() <= endDate: | |
|
|
694 |
|
|
|
|
692 | startDateTime= datetime.datetime.combine(startDate,startTime) | |
|
|
693 | endDateTime = datetime.datetime.combine(endDate,endTime) | |
|
|
694 | if startDateTime <= dt <= endDateTime: | |
|
|
695 | 695 |
|
|
|
696 | 696 | return False |
|
|
697 | 697 | |
| @@ -347,7 +347,7 class AMISRReader(ProcessingUnit): | |||
|
|
347 | 347 | |
|
|
348 | 348 | else: |
|
|
349 | 349 | #get the last file - 1 |
|
|
350 |
self.filenameList = [self.filenameList[- |
|
|
|
350 | self.filenameList = [self.filenameList[-2]] | |
|
|
351 | 351 | new_dirnameList = [] |
|
|
352 | 352 | for dirname in self.dirnameList: |
|
|
353 | 353 | junk = numpy.array([dirname in x for x in self.filenameList]) |
| @@ -421,7 +421,7 class AMISRReader(ProcessingUnit): | |||
|
|
421 | 421 | self.__selectDataForTimes(online=True) |
|
|
422 | 422 | filename = self.filenameList[0] |
|
|
423 | 423 | wait = 0 |
|
|
424 |
|
|
|
|
424 | self.__waitForNewFile=300 ## DEBUG: | |
|
|
425 | 425 | while self.__filename_online == filename: |
|
|
426 | 426 | print('waiting %d seconds to get a new file...'%(self.__waitForNewFile)) |
|
|
427 | 427 | if wait == 5: |
| @@ -641,13 +641,14 class AMISRReader(ProcessingUnit): | |||
|
|
641 | 641 | #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp |
|
|
642 | 642 | #print(self.dataOut.utctime) |
|
|
643 | 643 | self.dataOut.profileIndex = self.profileIndex |
|
|
644 | #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime) | |
|
|
644 | 645 | self.dataOut.flagNoData = False |
|
|
645 | 646 | # if indexprof == 0: |
|
|
646 | 647 | # print self.dataOut.utctime |
|
|
647 | 648 | |
|
|
648 | 649 | self.profileIndex += 1 |
|
|
649 | 650 | |
|
|
650 | return self.dataOut.data | |
|
|
651 | #return self.dataOut.data | |
|
|
651 | 652 | |
|
|
652 | 653 | |
|
|
653 | 654 | def run(self, **kwargs): |
| @@ -660,3 +661,4 class AMISRReader(ProcessingUnit): | |||
|
|
660 | 661 | self.isConfig = True |
|
|
661 | 662 | |
|
|
662 | 663 | self.getData() |
|
|
664 | return(self.dataOut.data) | |
| @@ -5,7 +5,6 from schainpy.model.data.jrodata import SpectraHeis | |||
|
|
5 | 5 | from schainpy.utils import log |
|
|
6 | 6 | |
|
|
7 | 7 | |
|
|
8 | ||
|
|
9 | 8 | class SpectraHeisProc(ProcessingUnit): |
|
|
10 | 9 | |
|
|
11 | 10 | def __init__(self):#, **kwargs): |
| This diff has been collapsed as it changes many lines, (559 lines changed) Show them Hide them | |||
| @@ -12,12 +12,15 import time | |||
|
|
12 | 12 | import itertools |
|
|
13 | 13 | |
|
|
14 | 14 | import numpy |
|
|
15 | import math | |
|
|
15 | 16 | |
|
|
16 | 17 | from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation |
|
|
17 | 18 | from schainpy.model.data.jrodata import Spectra |
|
|
18 | 19 | from schainpy.model.data.jrodata import hildebrand_sekhon |
|
|
19 | 20 | from schainpy.utils import log |
|
|
20 | 21 | |
|
|
22 | from scipy.optimize import curve_fit | |
|
|
23 | ||
|
|
21 | 24 | |
|
|
22 | 25 | class SpectraProc(ProcessingUnit): |
|
|
23 | 26 | |
| @@ -478,6 +481,562 class removeDC(Operation): | |||
|
|
478 | 481 | |
|
|
479 | 482 | return self.dataOut |
|
|
480 | 483 | |
|
|
484 | # import matplotlib.pyplot as plt | |
|
|
485 | ||
|
|
486 | def fit_func( x, a0, a1, a2): #, a3, a4, a5): | |
|
|
487 | z = (x - a1) / a2 | |
|
|
488 | y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2 | |
|
|
489 | return y | |
|
|
490 | class CleanRayleigh(Operation): | |
|
|
491 | ||
|
|
492 | def __init__(self): | |
|
|
493 | ||
|
|
494 | Operation.__init__(self) | |
|
|
495 | self.i=0 | |
|
|
496 | self.isConfig = False | |
|
|
497 | self.__dataReady = False | |
|
|
498 | self.__profIndex = 0 | |
|
|
499 | self.byTime = False | |
|
|
500 | self.byProfiles = False | |
|
|
501 | ||
|
|
502 | self.bloques = None | |
|
|
503 | self.bloque0 = None | |
|
|
504 | ||
|
|
505 | self.index = 0 | |
|
|
506 | ||
|
|
507 | self.buffer = 0 | |
|
|
508 | self.buffer2 = 0 | |
|
|
509 | self.buffer3 = 0 | |
|
|
510 | #self.min_hei = None | |
|
|
511 | #self.max_hei = None | |
|
|
512 | ||
|
|
513 | def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv): | |
|
|
514 | ||
|
|
515 | self.nChannels = dataOut.nChannels | |
|
|
516 | self.nProf = dataOut.nProfiles | |
|
|
517 | self.nPairs = dataOut.data_cspc.shape[0] | |
|
|
518 | self.pairsArray = numpy.array(dataOut.pairsList) | |
|
|
519 | self.spectra = dataOut.data_spc | |
|
|
520 | self.cspectra = dataOut.data_cspc | |
|
|
521 | self.heights = dataOut.heightList #alturas totales | |
|
|
522 | self.nHeights = len(self.heights) | |
|
|
523 | self.min_hei = min_hei | |
|
|
524 | self.max_hei = max_hei | |
|
|
525 | if (self.min_hei == None): | |
|
|
526 | self.min_hei = 0 | |
|
|
527 | if (self.max_hei == None): | |
|
|
528 | self.max_hei = dataOut.heightList[-1] | |
|
|
529 | self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero() | |
|
|
530 | self.heightsClean = self.heights[self.hval] #alturas filtradas | |
|
|
531 | self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas | |
|
|
532 | self.nHeightsClean = len(self.heightsClean) | |
|
|
533 | self.channels = dataOut.channelList | |
|
|
534 | self.nChan = len(self.channels) | |
|
|
535 | self.nIncohInt = dataOut.nIncohInt | |
|
|
536 | self.__initime = dataOut.utctime | |
|
|
537 | self.maxAltInd = self.hval[-1]+1 | |
|
|
538 | self.minAltInd = self.hval[0] | |
|
|
539 | ||
|
|
540 | self.crosspairs = dataOut.pairsList | |
|
|
541 | self.nPairs = len(self.crosspairs) | |
|
|
542 | self.normFactor = dataOut.normFactor | |
|
|
543 | self.nFFTPoints = dataOut.nFFTPoints | |
|
|
544 | self.ippSeconds = dataOut.ippSeconds | |
|
|
545 | self.currentTime = self.__initime | |
|
|
546 | self.pairsArray = numpy.array(dataOut.pairsList) | |
|
|
547 | self.factor_stdv = factor_stdv | |
|
|
548 | ||
|
|
549 | ||
|
|
550 | ||
|
|
551 | if n != None : | |
|
|
552 | self.byProfiles = True | |
|
|
553 | self.nIntProfiles = n | |
|
|
554 | else: | |
|
|
555 | self.__integrationtime = timeInterval | |
|
|
556 | ||
|
|
557 | self.__dataReady = False | |
|
|
558 | self.isConfig = True | |
|
|
559 | ||
|
|
560 | ||
|
|
561 | ||
|
|
562 | def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5): | |
|
|
563 | #print (dataOut.utctime) | |
|
|
564 | if not self.isConfig : | |
|
|
565 | #print("Setting config") | |
|
|
566 | self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv) | |
|
|
567 | #print("Config Done") | |
|
|
568 | tini=dataOut.utctime | |
|
|
569 | ||
|
|
570 | if self.byProfiles: | |
|
|
571 | if self.__profIndex == self.nIntProfiles: | |
|
|
572 | self.__dataReady = True | |
|
|
573 | else: | |
|
|
574 | if (tini - self.__initime) >= self.__integrationtime: | |
|
|
575 | #print(tini - self.__initime,self.__profIndex) | |
|
|
576 | self.__dataReady = True | |
|
|
577 | self.__initime = tini | |
|
|
578 | ||
|
|
579 | #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0): | |
|
|
580 | ||
|
|
581 | if self.__dataReady: | |
|
|
582 | print("Data ready",self.__profIndex) | |
|
|
583 | self.__profIndex = 0 | |
|
|
584 | jspc = self.buffer | |
|
|
585 | jcspc = self.buffer2 | |
|
|
586 | #jnoise = self.buffer3 | |
|
|
587 | self.buffer = dataOut.data_spc | |
|
|
588 | self.buffer2 = dataOut.data_cspc | |
|
|
589 | #self.buffer3 = dataOut.noise | |
|
|
590 | self.currentTime = dataOut.utctime | |
|
|
591 | if numpy.any(jspc) : | |
|
|
592 | #print( jspc.shape, jcspc.shape) | |
|
|
593 | jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights)) | |
|
|
594 | jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) | |
|
|
595 | self.__dataReady = False | |
|
|
596 | #print( jspc.shape, jcspc.shape) | |
|
|
597 | dataOut.flagNoData = False | |
|
|
598 | else: | |
|
|
599 | dataOut.flagNoData = True | |
|
|
600 | self.__dataReady = False | |
|
|
601 | return dataOut | |
|
|
602 | else: | |
|
|
603 | #print( len(self.buffer)) | |
|
|
604 | if numpy.any(self.buffer): | |
|
|
605 | self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0) | |
|
|
606 | self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) | |
|
|
607 | self.buffer3 += dataOut.data_dc | |
|
|
608 | else: | |
|
|
609 | self.buffer = dataOut.data_spc | |
|
|
610 | self.buffer2 = dataOut.data_cspc | |
|
|
611 | self.buffer3 = dataOut.data_dc | |
|
|
612 | #print self.index, self.fint | |
|
|
613 | #print self.buffer2.shape | |
|
|
614 | dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO | |
|
|
615 | self.__profIndex += 1 | |
|
|
616 | return dataOut ## NOTE: REV | |
|
|
617 | ||
|
|
618 | # if self.index == 0 and self.fint == 1 : | |
|
|
619 | # if jspc != None: | |
|
|
620 | # print len(jspc), jspc.shape | |
|
|
621 | # jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4)) | |
|
|
622 | # print jspc.shape | |
|
|
623 | # dataOut.flagNoData = True | |
|
|
624 | # return dataOut | |
|
|
625 | # if path != None: | |
|
|
626 | # sys.path.append(path) | |
|
|
627 | # self.library = importlib.import_module(file) | |
|
|
628 | # | |
|
|
629 | # To be inserted as a parameter | |
|
|
630 | ||
|
|
631 | #Set constants | |
|
|
632 | #constants = self.library.setConstants(dataOut) | |
|
|
633 | #dataOut.constants = constants | |
|
|
634 | ||
|
|
635 | #snrth= 20 | |
|
|
636 | ||
|
|
637 | #crosspairs = dataOut.groupList | |
|
|
638 | #noise = dataOut.noise | |
|
|
639 | #print( nProf,heights) | |
|
|
640 | #print( jspc.shape, jspc.shape[0]) | |
|
|
641 | #print noise | |
|
|
642 | #print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise) | |
|
|
643 | #jnoise = jnoise/N | |
|
|
644 | #noise = numpy.nansum(jnoise,axis=0)#/len(jnoise) | |
|
|
645 | #print( noise) | |
|
|
646 | #power = numpy.sum(spectra, axis=1) | |
|
|
647 | #print power[0,:] | |
|
|
648 | #print("CROSSPAIRS",crosspairs) | |
|
|
649 | #nPairs = len(crosspairs) | |
|
|
650 | #print(numpy.shape(dataOut.data_spc)) | |
|
|
651 | ||
|
|
652 | #absc = dataOut.abscissaList[:-1] | |
|
|
653 | ||
|
|
654 | #print absc.shape | |
|
|
655 | #nBlocks=149 | |
|
|
656 | #print('spectra', spectra.shape) | |
|
|
657 | #print('noise print', crosspairs) | |
|
|
658 | #print('spectra', spectra.shape) | |
|
|
659 | #print('cspectra', cspectra.shape) | |
|
|
660 | #print numpy.array(dataOut.data_pre[1]).shape | |
|
|
661 | #spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs) | |
|
|
662 | ||
|
|
663 | ||
|
|
664 | #index = tini.tm_hour*12+tini.tm_min/5 | |
|
|
665 | # jspc = jspc/self.nFFTPoints/self.normFactor | |
|
|
666 | # jcspc = jcspc/self.nFFTPoints/self.normFactor | |
|
|
667 | ||
|
|
668 | ||
|
|
669 | ||
|
|
670 | #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks) | |
|
|
671 | #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei) | |
|
|
672 | tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) | |
|
|
673 | #jspectra = tmp_spectra*len(jspc[:,0,0,0]) | |
|
|
674 | #jcspectra = tmp_cspectra*len(jspc[:,0,0,0]) | |
|
|
675 | ||
|
|
676 | dataOut.data_spc = tmp_spectra | |
|
|
677 | dataOut.data_cspc = tmp_cspectra | |
|
|
678 | dataOut.data_dc = self.buffer3 | |
|
|
679 | dataOut.nIncohInt *= self.nIntProfiles | |
|
|
680 | dataOut.utctime = self.currentTime #tiempo promediado | |
|
|
681 | print("Time: ",time.localtime(dataOut.utctime)) | |
|
|
682 | # dataOut.data_spc = sat_spectra | |
|
|
683 | # dataOut.data_cspc = sat_cspectra | |
|
|
684 | self.buffer = 0 | |
|
|
685 | self.buffer2 = 0 | |
|
|
686 | self.buffer3 = 0 | |
|
|
687 | ||
|
|
688 | return dataOut | |
|
|
689 | ||
|
|
690 | def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv): | |
|
|
691 | print("OP cleanRayleigh") | |
|
|
692 | #import matplotlib.pyplot as plt | |
|
|
693 | #for k in range(149): | |
|
|
694 | ||
|
|
695 | rfunc = cspectra.copy() #self.bloques | |
|
|
696 | val_spc = spectra*0.0 #self.bloque0*0.0 | |
|
|
697 | val_cspc = cspectra*0.0 #self.bloques*0.0 | |
|
|
698 | in_sat_spectra = spectra.copy() #self.bloque0 | |
|
|
699 | in_sat_cspectra = cspectra.copy() #self.bloques | |
|
|
700 | ||
|
|
701 | raxs = math.ceil(math.sqrt(self.nPairs)) | |
|
|
702 | caxs = math.ceil(self.nPairs/raxs) | |
|
|
703 | ||
|
|
704 | #print(self.hval) | |
|
|
705 | #print numpy.absolute(rfunc[:,0,0,14]) | |
|
|
706 | for ih in range(self.minAltInd,self.maxAltInd): | |
|
|
707 | for ifreq in range(self.nFFTPoints): | |
|
|
708 | # fig, axs = plt.subplots(raxs, caxs) | |
|
|
709 | # fig2, axs2 = plt.subplots(raxs, caxs) | |
|
|
710 | col_ax = 0 | |
|
|
711 | row_ax = 0 | |
|
|
712 | for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS | |
|
|
713 | #print("ii: ",ii) | |
|
|
714 | if (col_ax%caxs==0 and col_ax!=0): | |
|
|
715 | col_ax = 0 | |
|
|
716 | row_ax += 1 | |
|
|
717 | func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia? | |
|
|
718 | #print(func2clean.shape) | |
|
|
719 | val = (numpy.isfinite(func2clean)==True).nonzero() | |
|
|
720 | ||
|
|
721 | if len(val)>0: #limitador | |
|
|
722 | min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40) | |
|
|
723 | if min_val <= -40 : | |
|
|
724 | min_val = -40 | |
|
|
725 | max_val = numpy.around(numpy.amax(func2clean)+2) #< 200 | |
|
|
726 | if max_val >= 200 : | |
|
|
727 | max_val = 200 | |
|
|
728 | #print min_val, max_val | |
|
|
729 | step = 1 | |
|
|
730 | #print("Getting bins and the histogram") | |
|
|
731 | x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step | |
|
|
732 | y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step)) | |
|
|
733 | #print(len(y_dist),len(binstep[:-1])) | |
|
|
734 | #print(row_ax,col_ax, " ..") | |
|
|
735 | #print(self.pairsArray[ii][0],self.pairsArray[ii][1]) | |
|
|
736 | mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist) | |
|
|
737 | sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist)) | |
|
|
738 | parg = [numpy.amax(y_dist),mean,sigma] | |
|
|
739 | gauss_fit, covariance = None, None | |
|
|
740 | newY = None | |
|
|
741 | try : | |
|
|
742 | gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg) | |
|
|
743 | mode = gauss_fit[1] | |
|
|
744 | stdv = gauss_fit[2] | |
|
|
745 | #print(" FIT OK",gauss_fit) | |
|
|
746 | ''' | |
|
|
747 | newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2]) | |
|
|
748 | axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green') | |
|
|
749 | axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red') | |
|
|
750 | axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))''' | |
|
|
751 | except: | |
|
|
752 | mode = mean | |
|
|
753 | stdv = sigma | |
|
|
754 | #print("FIT FAIL") | |
|
|
755 | ||
|
|
756 | ||
|
|
757 | #print(mode,stdv) | |
|
|
758 | #Removing echoes greater than mode + 3*stdv | |
|
|
759 | #factor_stdv = 2 | |
|
|
760 | noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero() | |
|
|
761 | #noval tiene los indices que se van a remover | |
|
|
762 | #print("Pair ",ii," novals: ",len(noval[0])) | |
|
|
763 | if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N) | |
|
|
764 | novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero() | |
|
|
765 | #print(novall) | |
|
|
766 | #print(" ") | |
|
|
767 | cross_pairs = self.pairsArray[ii] | |
|
|
768 | #Getting coherent echoes which are removed. | |
|
|
769 | if len(novall[0]) > 0: | |
|
|
770 | ||
|
|
771 | val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1 | |
|
|
772 | val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1 | |
|
|
773 | val_cspc[novall[0],ii,ifreq,ih] = 1 | |
|
|
774 | #print("OUT NOVALL 1") | |
|
|
775 | #Removing coherent from ISR data | |
|
|
776 | ||
|
|
777 | #print(spectra[:,ii,ifreq,ih]) | |
|
|
778 | new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0]) | |
|
|
779 | mean_cspc = numpy.mean(new_a) | |
|
|
780 | new_b = numpy.delete(spectra[:,cross_pairs[0],ifreq,ih], noval[0]) | |
|
|
781 | mean_spc0 = numpy.mean(new_b) | |
|
|
782 | new_c = numpy.delete(spectra[:,cross_pairs[1],ifreq,ih], noval[0]) | |
|
|
783 | mean_spc1 = numpy.mean(new_c) | |
|
|
784 | spectra[noval,cross_pairs[0],ifreq,ih] = mean_spc0 | |
|
|
785 | spectra[noval,cross_pairs[1],ifreq,ih] = mean_spc1 | |
|
|
786 | cspectra[noval,ii,ifreq,ih] = mean_cspc | |
|
|
787 | ||
|
|
788 | ''' | |
|
|
789 | func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih])) | |
|
|
790 | y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step)) | |
|
|
791 | axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red') | |
|
|
792 | axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green') | |
|
|
793 | axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii])) | |
|
|
794 | ''' | |
|
|
795 | ||
|
|
796 | col_ax += 1 #contador de ploteo columnas | |
|
|
797 | ##print(col_ax) | |
|
|
798 | ''' | |
|
|
799 | title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km" | |
|
|
800 | title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED" | |
|
|
801 | fig.suptitle(title) | |
|
|
802 | fig2.suptitle(title2) | |
|
|
803 | plt.show()''' | |
|
|
804 | ||
|
|
805 | ''' channels = channels | |
|
|
806 | cross_pairs = cross_pairs | |
|
|
807 | #print("OUT NOVALL 2") | |
|
|
808 | ||
|
|
809 | vcross0 = (cross_pairs[0] == channels[ii]).nonzero() | |
|
|
810 | vcross1 = (cross_pairs[1] == channels[ii]).nonzero() | |
|
|
811 | vcross = numpy.concatenate((vcross0,vcross1),axis=None) | |
|
|
812 | #print('vcros =', vcross) | |
|
|
813 | ||
|
|
814 | #Getting coherent echoes which are removed. | |
|
|
815 | if len(novall) > 0: | |
|
|
816 | #val_spc[novall,ii,ifreq,ih] = 1 | |
|
|
817 | val_spc[ii,ifreq,ih,novall] = 1 | |
|
|
818 | if len(vcross) > 0: | |
|
|
819 | val_cspc[vcross,ifreq,ih,novall] = 1 | |
|
|
820 | ||
|
|
821 | #Removing coherent from ISR data. | |
|
|
822 | self.bloque0[ii,ifreq,ih,noval] = numpy.nan | |
|
|
823 | if len(vcross) > 0: | |
|
|
824 | self.bloques[vcross,ifreq,ih,noval] = numpy.nan | |
|
|
825 | ''' | |
|
|
826 | ||
|
|
827 | print("Getting average of the spectra and cross-spectra from incoherent echoes.") | |
|
|
828 | out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan | |
|
|
829 | out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan | |
|
|
830 | for ih in range(self.nHeights): | |
|
|
831 | for ifreq in range(self.nFFTPoints): | |
|
|
832 | for ich in range(self.nChan): | |
|
|
833 | tmp = spectra[:,ich,ifreq,ih] | |
|
|
834 | valid = (numpy.isfinite(tmp[:])==True).nonzero() | |
|
|
835 | # if ich == 0 and ifreq == 0 and ih == 17 : | |
|
|
836 | # print tmp | |
|
|
837 | # print valid | |
|
|
838 | # print len(valid[0]) | |
|
|
839 | #print('TMP',tmp) | |
|
|
840 | if len(valid[0]) >0 : | |
|
|
841 | out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) | |
|
|
842 | #for icr in range(nPairs): | |
|
|
843 | for icr in range(self.nPairs): | |
|
|
844 | tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih]) | |
|
|
845 | valid = (numpy.isfinite(tmp)==True).nonzero() | |
|
|
846 | if len(valid[0]) > 0: | |
|
|
847 | out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) | |
|
|
848 | ''' | |
|
|
849 | # print('##########################################################') | |
|
|
850 | print("Removing fake coherent echoes (at least 4 points around the point)") | |
|
|
851 | ||
|
|
852 | val_spectra = numpy.sum(val_spc,0) | |
|
|
853 | val_cspectra = numpy.sum(val_cspc,0) | |
|
|
854 | ||
|
|
855 | val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4) | |
|
|
856 | val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4) | |
|
|
857 | ||
|
|
858 | for i in range(nChan): | |
|
|
859 | for j in range(nProf): | |
|
|
860 | for k in range(nHeights): | |
|
|
861 | if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 : | |
|
|
862 | val_spc[:,i,j,k] = 0.0 | |
|
|
863 | for i in range(nPairs): | |
|
|
864 | for j in range(nProf): | |
|
|
865 | for k in range(nHeights): | |
|
|
866 | if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 : | |
|
|
867 | val_cspc[:,i,j,k] = 0.0 | |
|
|
868 | ||
|
|
869 | # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan)) | |
|
|
870 | # if numpy.isfinite(val_spectra)==str(True): | |
|
|
871 | # noval = (val_spectra<1).nonzero() | |
|
|
872 | # if len(noval) > 0: | |
|
|
873 | # val_spc[:,noval] = 0.0 | |
|
|
874 | # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights)) | |
|
|
875 | ||
|
|
876 | #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf)) | |
|
|
877 | #if numpy.isfinite(val_cspectra)==str(True): | |
|
|
878 | # noval = (val_cspectra<1).nonzero() | |
|
|
879 | # if len(noval) > 0: | |
|
|
880 | # val_cspc[:,noval] = 0.0 | |
|
|
881 | # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights)) | |
|
|
882 | tmp_sat_spectra = spectra.copy() | |
|
|
883 | tmp_sat_spectra = tmp_sat_spectra*numpy.nan | |
|
|
884 | tmp_sat_cspectra = cspectra.copy() | |
|
|
885 | tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan | |
|
|
886 | ''' | |
|
|
887 | # fig = plt.figure(figsize=(6,5)) | |
|
|
888 | # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 | |
|
|
889 | # ax = fig.add_axes([left, bottom, width, height]) | |
|
|
890 | # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:]))) | |
|
|
891 | # ax.clabel(cp, inline=True,fontsize=10) | |
|
|
892 | # plt.show() | |
|
|
893 | ''' | |
|
|
894 | val = (val_spc > 0).nonzero() | |
|
|
895 | if len(val[0]) > 0: | |
|
|
896 | tmp_sat_spectra[val] = in_sat_spectra[val] | |
|
|
897 | val = (val_cspc > 0).nonzero() | |
|
|
898 | if len(val[0]) > 0: | |
|
|
899 | tmp_sat_cspectra[val] = in_sat_cspectra[val] | |
|
|
900 | ||
|
|
901 | print("Getting average of the spectra and cross-spectra from incoherent echoes 2") | |
|
|
902 | sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float) | |
|
|
903 | sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex) | |
|
|
904 | for ih in range(nHeights): | |
|
|
905 | for ifreq in range(nProf): | |
|
|
906 | for ich in range(nChan): | |
|
|
907 | tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih]) | |
|
|
908 | valid = (numpy.isfinite(tmp)).nonzero() | |
|
|
909 | if len(valid[0]) > 0: | |
|
|
910 | sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) | |
|
|
911 | ||
|
|
912 | for icr in range(nPairs): | |
|
|
913 | tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih]) | |
|
|
914 | valid = (numpy.isfinite(tmp)).nonzero() | |
|
|
915 | if len(valid[0]) > 0: | |
|
|
916 | sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) | |
|
|
917 | ''' | |
|
|
918 | #self.__dataReady= True | |
|
|
919 | #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra | |
|
|
920 | #if not self.__dataReady: | |
|
|
921 | #return None, None | |
|
|
922 | #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra | |
|
|
923 | return out_spectra, out_cspectra | |
|
|
924 | ||
|
|
925 | def REM_ISOLATED_POINTS(self,array,rth): | |
|
|
926 | # import matplotlib.pyplot as plt | |
|
|
927 | if rth == None : | |
|
|
928 | rth = 4 | |
|
|
929 | print("REM ISO") | |
|
|
930 | num_prof = len(array[0,:,0]) | |
|
|
931 | num_hei = len(array[0,0,:]) | |
|
|
932 | n2d = len(array[:,0,0]) | |
|
|
933 | ||
|
|
934 | for ii in range(n2d) : | |
|
|
935 | #print ii,n2d | |
|
|
936 | tmp = array[ii,:,:] | |
|
|
937 | #print tmp.shape, array[ii,101,:],array[ii,102,:] | |
|
|
938 | ||
|
|
939 | # fig = plt.figure(figsize=(6,5)) | |
|
|
940 | # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 | |
|
|
941 | # ax = fig.add_axes([left, bottom, width, height]) | |
|
|
942 | # x = range(num_prof) | |
|
|
943 | # y = range(num_hei) | |
|
|
944 | # cp = ax.contour(y,x,tmp) | |
|
|
945 | # ax.clabel(cp, inline=True,fontsize=10) | |
|
|
946 | # plt.show() | |
|
|
947 | ||
|
|
948 | #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs) | |
|
|
949 | tmp = numpy.reshape(tmp,num_prof*num_hei) | |
|
|
950 | indxs1 = (numpy.isfinite(tmp)==True).nonzero() | |
|
|
951 | indxs2 = (tmp > 0).nonzero() | |
|
|
952 | ||
|
|
953 | indxs1 = (indxs1[0]) | |
|
|
954 | indxs2 = indxs2[0] | |
|
|
955 | #indxs1 = numpy.array(indxs1[0]) | |
|
|
956 | #indxs2 = numpy.array(indxs2[0]) | |
|
|
957 | indxs = None | |
|
|
958 | #print indxs1 , indxs2 | |
|
|
959 | for iv in range(len(indxs2)): | |
|
|
960 | indv = numpy.array((indxs1 == indxs2[iv]).nonzero()) | |
|
|
961 | #print len(indxs2), indv | |
|
|
962 | if len(indv[0]) > 0 : | |
|
|
963 | indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None) | |
|
|
964 | # print indxs | |
|
|
965 | indxs = indxs[1:] | |
|
|
966 | #print(indxs, len(indxs)) | |
|
|
967 | if len(indxs) < 4 : | |
|
|
968 | array[ii,:,:] = 0. | |
|
|
969 | return | |
|
|
970 | ||
|
|
971 | xpos = numpy.mod(indxs ,num_hei) | |
|
|
972 | ypos = (indxs / num_hei) | |
|
|
973 | sx = numpy.argsort(xpos) # Ordering respect to "x" (time) | |
|
|
974 | #print sx | |
|
|
975 | xpos = xpos[sx] | |
|
|
976 | ypos = ypos[sx] | |
|
|
977 | ||
|
|
978 | # *********************************** Cleaning isolated points ********************************** | |
|
|
979 | ic = 0 | |
|
|
980 | while True : | |
|
|
981 | r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2))) | |
|
|
982 | #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh) | |
|
|
983 | #plt.plot(r) | |
|
|
984 | #plt.show() | |
|
|
985 | no_coh1 = (numpy.isfinite(r)==True).nonzero() | |
|
|
986 | no_coh2 = (r <= rth).nonzero() | |
|
|
987 | #print r, no_coh1, no_coh2 | |
|
|
988 | no_coh1 = numpy.array(no_coh1[0]) | |
|
|
989 | no_coh2 = numpy.array(no_coh2[0]) | |
|
|
990 | no_coh = None | |
|
|
991 | #print valid1 , valid2 | |
|
|
992 | for iv in range(len(no_coh2)): | |
|
|
993 | indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero()) | |
|
|
994 | if len(indv[0]) > 0 : | |
|
|
995 | no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None) | |
|
|
996 | no_coh = no_coh[1:] | |
|
|
997 | #print len(no_coh), no_coh | |
|
|
998 | if len(no_coh) < 4 : | |
|
|
999 | #print xpos[ic], ypos[ic], ic | |
|
|
1000 | # plt.plot(r) | |
|
|
1001 | # plt.show() | |
|
|
1002 | xpos[ic] = numpy.nan | |
|
|
1003 | ypos[ic] = numpy.nan | |
|
|
1004 | ||
|
|
1005 | ic = ic + 1 | |
|
|
1006 | if (ic == len(indxs)) : | |
|
|
1007 | break | |
|
|
1008 | #print( xpos, ypos) | |
|
|
1009 | ||
|
|
1010 | indxs = (numpy.isfinite(list(xpos))==True).nonzero() | |
|
|
1011 | #print indxs[0] | |
|
|
1012 | if len(indxs[0]) < 4 : | |
|
|
1013 | array[ii,:,:] = 0. | |
|
|
1014 | return | |
|
|
1015 | ||
|
|
1016 | xpos = xpos[indxs[0]] | |
|
|
1017 | ypos = ypos[indxs[0]] | |
|
|
1018 | for i in range(0,len(ypos)): | |
|
|
1019 | ypos[i]=int(ypos[i]) | |
|
|
1020 | junk = tmp | |
|
|
1021 | tmp = junk*0.0 | |
|
|
1022 | ||
|
|
1023 | tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))] | |
|
|
1024 | array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei)) | |
|
|
1025 | ||
|
|
1026 | #print array.shape | |
|
|
1027 | #tmp = numpy.reshape(tmp,(num_prof,num_hei)) | |
|
|
1028 | #print tmp.shape | |
|
|
1029 | ||
|
|
1030 | # fig = plt.figure(figsize=(6,5)) | |
|
|
1031 | # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 | |
|
|
1032 | # ax = fig.add_axes([left, bottom, width, height]) | |
|
|
1033 | # x = range(num_prof) | |
|
|
1034 | # y = range(num_hei) | |
|
|
1035 | # cp = ax.contour(y,x,array[ii,:,:]) | |
|
|
1036 | # ax.clabel(cp, inline=True,fontsize=10) | |
|
|
1037 | # plt.show() | |
|
|
1038 | return array | |
|
|
1039 | ||
|
|
481 | 1040 | class removeInterference(Operation): |
|
|
482 | 1041 | |
|
|
483 | 1042 | def removeInterference2(self): |
General Comments 0
You need to be logged in to leave comments.
Login now