| @@ -3158,10 +3158,59 class SpectralFitting(Operation): | |||
|
|
3158 | 3158 | self.isConfig = False |
|
|
3159 | 3159 | |
|
|
3160 | 3160 | |
|
|
3161 | def setup(self,nChan,nProf,nHei,nBlocks): | |
|
|
3161 | def setup(self,dataOut,groupList,path,file,filec): | |
|
|
3162 | 3162 | self.__dataReady = False |
|
|
3163 | self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex) | |
|
|
3164 | self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks]) | |
|
|
3163 | # new | |
|
|
3164 | self.nChannels = dataOut.nChannels | |
|
|
3165 | self.channels = dataOut.channelList | |
|
|
3166 | self.nHeights = dataOut.heightList.size | |
|
|
3167 | self.heights = dataOut.heightList | |
|
|
3168 | self.nProf = dataOut.nProfiles | |
|
|
3169 | self.nIncohInt = dataOut.nIncohInt | |
|
|
3170 | self.absc = dataOut.abscissaList[:-1] | |
|
|
3171 | ||
|
|
3172 | ||
|
|
3173 | #To be inserted as a parameter | |
|
|
3174 | try: | |
|
|
3175 | self.groupArray = numpy.array(groupList)#groupArray = numpy.array([[0,1],[2,3]]) | |
|
|
3176 | except: | |
|
|
3177 | print("Please insert groupList. Example (0,1),(2,3) format multilist") | |
|
|
3178 | dataOut.groupList = self.groupArray | |
|
|
3179 | self.crosspairs = dataOut.groupList | |
|
|
3180 | self.nPairs = len(self.crosspairs) | |
|
|
3181 | self.nGroups = self.groupArray.shape[0] | |
|
|
3182 | ||
|
|
3183 | #List of possible combinations | |
|
|
3184 | ||
|
|
3185 | self.listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2) | |
|
|
3186 | self.indCross = numpy.zeros(len(list(self.listComb)), dtype = 'int') | |
|
|
3187 | ||
|
|
3188 | #Parameters Array | |
|
|
3189 | dataOut.data_param = None | |
|
|
3190 | dataOut.data_paramC = None | |
|
|
3191 | dataOut.clean_num_aver = None | |
|
|
3192 | dataOut.coh_num_aver = None | |
|
|
3193 | dataOut.tmp_spectra_i = None | |
|
|
3194 | dataOut.tmp_cspectra_i = None | |
|
|
3195 | dataOut.tmp_spectra_c = None | |
|
|
3196 | dataOut.tmp_cspectra_c = None | |
|
|
3197 | dataOut.index = None | |
|
|
3198 | ||
|
|
3199 | if path != None: | |
|
|
3200 | sys.path.append(path) | |
|
|
3201 | self.library = importlib.import_module(file) | |
|
|
3202 | if filec != None: | |
|
|
3203 | self.weightf = importlib.import_module(filec) | |
|
|
3204 | ||
|
|
3205 | #Set constants | |
|
|
3206 | self.constants = self.library.setConstants(dataOut) | |
|
|
3207 | dataOut.constants = self.constants | |
|
|
3208 | self.M = dataOut.normFactor | |
|
|
3209 | self.N = dataOut.nFFTPoints | |
|
|
3210 | self.ippSeconds = dataOut.ippSeconds | |
|
|
3211 | self.K = dataOut.nIncohInt | |
|
|
3212 | self.pairsArray = numpy.array(dataOut.pairsList) | |
|
|
3213 | self.snrth = 20 | |
|
|
3165 | 3214 | |
|
|
3166 | 3215 | def __calculateMoments(self,oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None): |
|
|
3167 | 3216 | |
| @@ -3343,7 +3392,6 class SpectralFitting(Operation): | |||
|
|
3343 | 3392 | incoh_aver[pair[1],incoh_echoes]=1 |
|
|
3344 | 3393 | return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver |
|
|
3345 | 3394 | |
|
|
3346 | ||
|
|
3347 | 3395 | def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index): |
|
|
3348 | 3396 | |
|
|
3349 | 3397 | nProf = dataOut.nProfiles |
| @@ -3520,6 +3568,7 class SpectralFitting(Operation): | |||
|
|
3520 | 3568 | if len(valid[0]) > 0: |
|
|
3521 | 3569 | sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) |
|
|
3522 | 3570 | return out_spectra, out_cspectra,sat_spectra,sat_cspectra |
|
|
3571 | ||
|
|
3523 | 3572 | def REM_ISOLATED_POINTS(self,array,rth): |
|
|
3524 | 3573 | if rth == None : rth = 4 |
|
|
3525 | 3574 | num_prof = len(array[0,:,0]) |
| @@ -3603,37 +3652,39 class SpectralFitting(Operation): | |||
|
|
3603 | 3652 | smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp) |
|
|
3604 | 3653 | return [fmom,numpy.sqrt(smom)] |
|
|
3605 | 3654 | |
|
|
3606 | ||
|
|
3607 | ||
|
|
3608 | ||
|
|
3609 | ||
|
|
3610 | 3655 | def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None): |
|
|
3611 | if not numpy.any(proc): | |
|
|
3612 | nChannels = dataOut.nChannels | |
|
|
3613 | nHeights= dataOut.heightList.size | |
|
|
3614 | nProf = dataOut.nProfiles | |
|
|
3615 |
|
|
|
|
3616 |
|
|
|
|
3617 | tini=time.localtime(dataOut.utctime) | |
|
|
3656 | if not self.isConfig: | |
|
|
3657 | self.setup(dataOut = dataOut,groupList=groupList,path=path,file=file,filec=filec) | |
|
|
3658 | self.isConfig = True | |
|
|
3659 | ||
|
|
3660 | if not numpy.any(proc): | |
|
|
3661 | if numpy.any(taver): | |
|
|
3662 | taver = int(taver) | |
|
|
3663 | else : | |
|
|
3664 | taver = 5 | |
|
|
3665 | tini = time.localtime(dataOut.utctime) | |
|
|
3618 | 3666 | if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0): |
|
|
3619 | self.index = 0 | |
|
|
3620 | jspc = self.buffer | |
|
|
3621 | jcspc = self.buffer2 | |
|
|
3622 | jnoise = self.buffer3 | |
|
|
3623 | self.buffer = dataOut.data_spc | |
|
|
3667 | self.index = 0 | |
|
|
3668 | jspc = self.buffer | |
|
|
3669 | jcspc = self.buffer2 | |
|
|
3670 | jnoise = self.buffer3 | |
|
|
3671 | self.buffer = dataOut.data_spc | |
|
|
3624 | 3672 | self.buffer2 = dataOut.data_cspc |
|
|
3625 | 3673 | self.buffer3 = dataOut.noise |
|
|
3626 | self.fint = 1 | |
|
|
3674 | self.fint = 1 | |
|
|
3627 | 3675 | if numpy.any(jspc) : |
|
|
3628 | jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights)) | |
|
|
3629 | jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights)) | |
|
|
3630 | jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels)) | |
|
|
3676 | jspc = numpy.reshape(jspc ,(int(len(jspc) / self.nChannels) , self.nChannels ,self.nProf,self.nHeights )) | |
|
|
3677 | jcspc = numpy.reshape(jcspc ,(int(len(jcspc) /int(self.nChannels/2)),int(self.nChannels/2),self.nProf,self.nHeights )) | |
|
|
3678 | jnoise = numpy.reshape(jnoise,(int(len(jnoise)/ self.nChannels) , self.nChannels)) | |
|
|
3631 | 3679 | else: |
|
|
3632 | 3680 | dataOut.flagNoData = True |
|
|
3633 | 3681 | return dataOut |
|
|
3634 | 3682 | else : |
|
|
3635 |
if (tini.tm_min % taver) == 0 : |
|
|
|
3636 |
|
|
|
|
3683 | if (tini.tm_min % taver) == 0 : | |
|
|
3684 | self.fint = 1 | |
|
|
3685 | else : | |
|
|
3686 | self.fint = 0 | |
|
|
3687 | ||
|
|
3637 | 3688 | self.index += 1 |
|
|
3638 | 3689 | if numpy.any(self.buffer): |
|
|
3639 | 3690 | self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0) |
| @@ -3645,85 +3696,36 class SpectralFitting(Operation): | |||
|
|
3645 | 3696 | self.buffer3 = dataOut.noise |
|
|
3646 | 3697 | dataOut.flagNoData = True |
|
|
3647 | 3698 | return dataOut |
|
|
3648 | if path != None: | |
|
|
3649 | sys.path.append(path) | |
|
|
3650 | self.library = importlib.import_module(file) | |
|
|
3651 | if filec != None: | |
|
|
3652 | self.weightf = importlib.import_module(filec) | |
|
|
3653 | ||
|
|
3654 | #To be inserted as a parameter | |
|
|
3655 | groupArray = numpy.array(groupList) | |
|
|
3656 | #groupArray = numpy.array([[0,1],[2,3]]) | |
|
|
3657 | dataOut.groupList = groupArray | |
|
|
3658 | nGroups = groupArray.shape[0] | |
|
|
3659 | nChannels = dataOut.nChannels | |
|
|
3660 | nHeights = dataOut.heightList.size | |
|
|
3661 | ||
|
|
3662 | #Parameters Array | |
|
|
3663 | dataOut.data_param = None | |
|
|
3664 | dataOut.data_paramC = None | |
|
|
3665 | dataOut.clean_num_aver = None | |
|
|
3666 | dataOut.coh_num_aver = None | |
|
|
3667 | dataOut.tmp_spectra_i = None | |
|
|
3668 | dataOut.tmp_cspectra_i = None | |
|
|
3669 | dataOut.tmp_spectra_c = None | |
|
|
3670 | dataOut.tmp_cspectra_c = None | |
|
|
3671 | dataOut.index = None | |
|
|
3672 | ||
|
|
3673 | #Set constants | |
|
|
3674 | constants = self.library.setConstants(dataOut) | |
|
|
3675 | dataOut.constants = constants | |
|
|
3676 | M = dataOut.normFactor | |
|
|
3677 | N = dataOut.nFFTPoints | |
|
|
3678 | ippSeconds = dataOut.ippSeconds | |
|
|
3679 | K = dataOut.nIncohInt | |
|
|
3680 | pairsArray = numpy.array(dataOut.pairsList) | |
|
|
3681 | snrth= 20 | |
|
|
3682 | spectra = dataOut.data_spc | |
|
|
3683 | cspectra = dataOut.data_cspc | |
|
|
3684 | nProf = dataOut.nProfiles | |
|
|
3685 | heights = dataOut.heightList | |
|
|
3686 | nHei = len(heights) | |
|
|
3687 | channels = dataOut.channelList | |
|
|
3688 | nChan = len(channels) | |
|
|
3689 | nIncohInt = dataOut.nIncohInt | |
|
|
3690 | crosspairs = dataOut.groupList | |
|
|
3691 | noise = dataOut.noise | |
|
|
3692 | jnoise = jnoise/N | |
|
|
3693 | noise = numpy.nansum(jnoise,axis=0)#/len(jnoise) | |
|
|
3694 | power = numpy.sum(spectra, axis=1) | |
|
|
3695 | nPairs = len(crosspairs) | |
|
|
3696 | absc = dataOut.abscissaList[:-1] | |
|
|
3697 | 3699 | |
|
|
3698 | if not self.isConfig: | |
|
|
3699 | self.isConfig = True | |
|
|
3700 | jnoise = jnoise/self.N# creo que falta dividirlo entre N | |
|
|
3701 | noise = numpy.nansum(jnoise,axis=0)#/len(jnoise) | |
|
|
3702 | index = tini.tm_hour*12+tini.tm_min/taver | |
|
|
3703 | dataOut.index = index | |
|
|
3704 | jspc = jspc/self.N/self.N | |
|
|
3705 | jcspc = jcspc/self.N/self.N | |
|
|
3700 | 3706 | |
|
|
3701 | index = tini.tm_hour*12+tini.tm_min/taver | |
|
|
3702 | dataOut.index= index | |
|
|
3703 | jspc = jspc/N/N | |
|
|
3704 | jcspc = jcspc/N/N | |
|
|
3705 | 3707 | tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2) |
|
|
3706 | jspectra = tmp_spectra*len(jspc[:,0,0,0]) | |
|
|
3707 | jcspectra = tmp_cspectra*len(jspc[:,0,0,0]) | |
|
|
3708 | my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th) | |
|
|
3709 | clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index) | |
|
|
3710 | dataOut.data_spc = incoh_spectra | |
|
|
3711 | dataOut.data_cspc = incoh_cspectra | |
|
|
3712 | clean_num_aver = incoh_aver*len(jspc[:,0,0,0]) | |
|
|
3713 | coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0]) | |
|
|
3714 | dataOut.clean_num_aver = clean_num_aver | |
|
|
3715 | dataOut.coh_num_aver = coh_num_aver | |
|
|
3716 | dataOut.tmp_spectra_i = incoh_spectra | |
|
|
3717 | dataOut.tmp_cspectra_i = incoh_cspectra | |
|
|
3718 | dataOut.tmp_spectra_c = clean_coh_spectra | |
|
|
3719 | dataOut.tmp_cspectra_c = clean_coh_cspectra | |
|
|
3708 | jspectra = tmp_spectra * len(jspc[:,0,0,0]) | |
|
|
3709 | jcspectra = tmp_cspectra * len(jspc[:,0,0,0]) | |
|
|
3710 | ||
|
|
3711 | my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, self.snrth,coh_th, hei_th) | |
|
|
3712 | clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(self.snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index) | |
|
|
3713 | ||
|
|
3714 | dataOut.data_spc = incoh_spectra | |
|
|
3715 | dataOut.data_cspc = incoh_cspectra | |
|
|
3716 | clean_num_aver = incoh_aver * len(jspc[:,0,0,0]) | |
|
|
3717 | coh_num_aver = clean_coh_aver* len(jspc[:,0,0,0]) | |
|
|
3718 | dataOut.clean_num_aver = clean_num_aver | |
|
|
3719 | dataOut.coh_num_aver = coh_num_aver | |
|
|
3720 | ||
|
|
3720 | 3721 | #List of possible combinations |
|
|
3721 | listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) | |
|
|
3722 | indCross = numpy.zeros(len(list(listComb)), dtype = 'int') | |
|
|
3723 | if getSNR: | |
|
|
3724 | listChannels = groupArray.reshape((groupArray.size)) | |
|
|
3725 | listChannels.sort() | |
|
|
3726 | dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) | |
|
|
3722 | #listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) | |
|
|
3723 | #indCross = numpy.zeros(len(list(listComb)), dtype = 'int') | |
|
|
3724 | #if getSNR: | |
|
|
3725 | # listChannels = groupArray.reshape((groupArray.size)) | |
|
|
3726 | # listChannels.sort() | |
|
|
3727 | # print("AQUI ESTOY") | |
|
|
3728 | # dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) | |
|
|
3727 | 3729 | else: |
|
|
3728 | 3730 | clean_num_aver = dataOut.clean_num_aver |
|
|
3729 | 3731 | coh_num_aver = dataOut.coh_num_aver |
| @@ -3762,51 +3764,49 class SpectralFitting(Operation): | |||
|
|
3762 | 3764 | constants['M'] = M |
|
|
3763 | 3765 | dataOut.constants = constants |
|
|
3764 | 3766 | |
|
|
3765 | groupArray = numpy.array(groupList) | |
|
|
3766 | dataOut.groupList = groupArray | |
|
|
3767 | nGroups = groupArray.shape[0] | |
|
|
3768 | 3767 | #List of possible combinations |
|
|
3769 | listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) | |
|
|
3768 | listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2) | |
|
|
3770 | 3769 | indCross = numpy.zeros(len(list(listComb)), dtype = 'int') |
|
|
3771 | 3770 | if dataOut.data_paramC is None: |
|
|
3772 |
|
|
|
|
3773 | for i in range(nGroups): | |
|
|
3774 | coord = groupArray[i,:] | |
|
|
3771 | dataOut.data_paramC = numpy.zeros((self.nGroups*4, self.nHeights ,2))*numpy.nan | |
|
|
3772 | for i in range(self.nGroups): | |
|
|
3773 | coord = self.groupArray[i,:] | |
|
|
3775 | 3774 | #Input data array |
|
|
3776 | data = dataOut.data_spc[coord,:,:]/(M*N) | |
|
|
3777 | data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) | |
|
|
3775 | data = dataOut.data_spc[coord,:,:]/(self.M*self.N) | |
|
|
3776 | data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) | |
|
|
3778 | 3777 | |
|
|
3779 | 3778 | #Cross Spectra data array for Covariance Matrixes |
|
|
3780 | 3779 | ind = 0 |
|
|
3781 | 3780 | for pairs in listComb: |
|
|
3782 | pairsSel = numpy.array([coord[x],coord[y]]) | |
|
|
3783 | indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) | |
|
|
3784 | ind += 1 | |
|
|
3785 | dataCross = dataOut.data_cspc[indCross,:,:]/(M*N) | |
|
|
3781 | pairsSel = numpy.array([coord[x],coord[y]]) | |
|
|
3782 | indCross[ind] = int(numpy.where(numpy.all(self.pairsArray == pairsSel, axis = 1))[0][0]) | |
|
|
3783 | ind += 1 | |
|
|
3784 | dataCross = dataOut.data_cspc[indCross,:,:]/(self.M*self.N) | |
|
|
3786 | 3785 | dataCross = dataCross**2 |
|
|
3787 | nhei = nHeights | |
|
|
3788 | poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:] | |
|
|
3786 | nhei = self.nHeights | |
|
|
3787 | poweri = numpy.sum(dataOut.data_spc[:,1:self.nProf-0,:],axis=1)/clean_num_aver[:,:] | |
|
|
3788 | ||
|
|
3789 | 3789 | if i == 0 : my_noises = numpy.zeros(4,dtype=float) |
|
|
3790 | n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1) | |
|
|
3791 | n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1) | |
|
|
3790 | n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(self.nProf-1) | |
|
|
3791 | n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(self.nProf-1) | |
|
|
3792 | 3792 | n0 = n0i |
|
|
3793 | 3793 | n1= n1i |
|
|
3794 | 3794 | my_noises[2*i+0] = n0 |
|
|
3795 | 3795 | my_noises[2*i+1] = n1 |
|
|
3796 | 3796 | snrth = -15.0 # -4 -16 -25 |
|
|
3797 | 3797 | snrth = 10**(snrth/10.0) |
|
|
3798 | jvelr = numpy.zeros(nHeights, dtype = 'float') | |
|
|
3798 | jvelr = numpy.zeros(self.nHeights, dtype = 'float') | |
|
|
3799 | 3799 | hvalid = [0] |
|
|
3800 | coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:]) | |
|
|
3800 | coh2 = abs(dataOut.data_cspc[i,1:self.nProf,:])**2/(dataOut.data_spc[0+i*2,1:self.nProf-0,:]*dataOut.data_spc[1+i*2,1:self.nProf-0,:]) | |
|
|
3801 | 3801 | |
|
|
3802 | for h in range(nHeights): | |
|
|
3802 | for h in range(self.nHeights): | |
|
|
3803 | 3803 | smooth = clean_num_aver[i+1,h] |
|
|
3804 | signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth | |
|
|
3805 | signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth | |
|
|
3804 | signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth | |
|
|
3805 | signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth | |
|
|
3806 | 3806 | signal0 = signalpn0-n0 |
|
|
3807 | 3807 | signal1 = signalpn1-n1 |
|
|
3808 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3809 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
3808 | snr0 = numpy.sum(signal0/n0)/(self.nProf-1) | |
|
|
3809 | snr1 = numpy.sum(signal1/n1)/(self.nProf-1) | |
|
|
3810 | 3810 | gamma = coh2[:,h] |
|
|
3811 | 3811 | indxs = (numpy.isfinite(list(gamma))==True).nonzero() |
|
|
3812 | 3812 | if len(indxs) >0: |
| @@ -3817,22 +3817,22 class SpectralFitting(Operation): | |||
|
|
3817 | 3817 | else: |
|
|
3818 | 3818 | maxp0 = numpy.argmax(signal0) |
|
|
3819 | 3819 | maxp1 = numpy.argmax(signal1) |
|
|
3820 | jvelr[h] = (absc[maxp0]+absc[maxp1])/2. | |
|
|
3821 | else: jvelr[h] = absc[0] | |
|
|
3820 | jvelr[h] = (self.absc[maxp0]+self.absc[maxp1])/2. | |
|
|
3821 | else: jvelr[h] = self.absc[0] | |
|
|
3822 | 3822 | if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None) |
|
|
3823 | 3823 | #print(maxp0,absc[maxp0],snr0,jvelr[h]) |
|
|
3824 | 3824 | |
|
|
3825 | 3825 | if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1 |
|
|
3826 | 3826 | else: fd0 = numpy.nan |
|
|
3827 | for h in range(nHeights): | |
|
|
3827 | for h in range(self.nHeights): | |
|
|
3828 | 3828 | d = data[:,h] |
|
|
3829 | 3829 | smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:] |
|
|
3830 | signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth | |
|
|
3831 | signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth | |
|
|
3830 | signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth | |
|
|
3831 | signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth | |
|
|
3832 | 3832 | signal0 = signalpn0-n0 |
|
|
3833 | 3833 | signal1 = signalpn1-n1 |
|
|
3834 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3835 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
3834 | snr0 = numpy.sum(signal0/n0)/(self.nProf-1) | |
|
|
3835 | snr1 = numpy.sum(signal1/n1)/(self.nProf-1) | |
|
|
3836 | 3836 | if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 : |
|
|
3837 | 3837 | #Covariance Matrix |
|
|
3838 | 3838 | D = numpy.diag(d**2) |
| @@ -3867,7 +3867,7 class SpectralFitting(Operation): | |||
|
|
3867 | 3867 | if (h>6)and(error1[3]<25): |
|
|
3868 | 3868 | p0 = dataOut.data_param[i,:,h-1].copy() |
|
|
3869 | 3869 | else: |
|
|
3870 | p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i) | |
|
|
3870 | p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, self.constants))# sin el i(data_spc, constants, i) | |
|
|
3871 | 3871 | p0[3] = fd0 |
|
|
3872 | 3872 | |
|
|
3873 | 3873 | if filec != None: |
| @@ -3875,12 +3875,12 class SpectralFitting(Operation): | |||
|
|
3875 | 3875 | |
|
|
3876 | 3876 | try: |
|
|
3877 | 3877 | #Least Squares |
|
|
3878 | minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) | |
|
|
3878 | minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,self.constants),full_output=True) | |
|
|
3879 | 3879 | #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) |
|
|
3880 | 3880 | #Chi square error |
|
|
3881 | error0 = numpy.sum(infodict['fvec']**2)/(2*N) | |
|
|
3881 | error0 = numpy.sum(infodict['fvec']**2)/(2*self.N) | |
|
|
3882 | 3882 | #Error with Jacobian |
|
|
3883 | error1 = self.library.errorFunction(minp,constants,LT) | |
|
|
3883 | error1 = self.library.errorFunction(minp,self.constants,LT) | |
|
|
3884 | 3884 | |
|
|
3885 | 3885 | except: |
|
|
3886 | 3886 | minp = p0*numpy.nan |
| @@ -3888,29 +3888,29 class SpectralFitting(Operation): | |||
|
|
3888 | 3888 | error1 = p0*numpy.nan |
|
|
3889 | 3889 | else : |
|
|
3890 | 3890 | data_spc = dataOut.data_spc[coord,:,h] |
|
|
3891 | p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants)) | |
|
|
3891 | p0 = numpy.array(self.library.initialValuesFunction(data_spc, self.constants)) | |
|
|
3892 | 3892 | minp = p0*numpy.nan |
|
|
3893 | 3893 | error0 = numpy.nan |
|
|
3894 | 3894 | error1 = p0*numpy.nan |
|
|
3895 | 3895 | if dataOut.data_param is None: |
|
|
3896 | dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan | |
|
|
3897 | dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan | |
|
|
3896 | dataOut.data_param = numpy.zeros((self.nGroups, p0.size, self.nHeights ))*numpy.nan | |
|
|
3897 | dataOut.data_error = numpy.zeros((self.nGroups, p0.size + 1, self.nHeights ))*numpy.nan | |
|
|
3898 | 3898 | dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) |
|
|
3899 | 3899 | dataOut.data_param[i,:,h] = minp |
|
|
3900 | 3900 | |
|
|
3901 | for ht in range(nHeights-1) : | |
|
|
3901 | for ht in range(self.nHeights-1) : | |
|
|
3902 | 3902 | smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam] |
|
|
3903 | 3903 | dataOut.data_paramC[4*i,ht,1] = smooth |
|
|
3904 | signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra | |
|
|
3905 | signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth | |
|
|
3904 | signalpn0 = (clean_coh_spectra[i*2 ,1:(self.nProf-0),ht])/smooth #coh_spectra | |
|
|
3905 | signalpn1 = (clean_coh_spectra[i*2+1,1:(self.nProf-0),ht])/smooth | |
|
|
3906 | 3906 | val0 = (signalpn0 > 0).nonzero() |
|
|
3907 | 3907 | val0 = val0[0] |
|
|
3908 | if len(val0) == 0 : val0_npoints = nProf | |
|
|
3908 | if len(val0) == 0 : val0_npoints = self.nProf | |
|
|
3909 | 3909 | else : val0_npoints = len(val0) |
|
|
3910 | 3910 | |
|
|
3911 | 3911 | val1 = (signalpn1 > 0).nonzero() |
|
|
3912 | 3912 | val1 = val1[0] |
|
|
3913 | if len(val1) == 0 : val1_npoints = nProf | |
|
|
3913 | if len(val1) == 0 : val1_npoints = self.nProf | |
|
|
3914 | 3914 | else : val1_npoints = len(val1) |
|
|
3915 | 3915 | |
|
|
3916 | 3916 | dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0 |
| @@ -3924,26 +3924,32 class SpectralFitting(Operation): | |||
|
|
3924 | 3924 | vali = (signal1 < 0).nonzero() |
|
|
3925 | 3925 | vali = vali[0] |
|
|
3926 | 3926 | if len(vali) > 0 : signal1[vali] = 0 |
|
|
3927 | snr0 = numpy.sum(signal0/n0)/(nProf-1) | |
|
|
3928 | snr1 = numpy.sum(signal1/n1)/(nProf-1) | |
|
|
3929 | doppler = absc[1:] | |
|
|
3927 | snr0 = numpy.sum(signal0/n0)/(self.nProf-1) | |
|
|
3928 | snr1 = numpy.sum(signal1/n1)/(self.nProf-1) | |
|
|
3929 | doppler = self.absc[1:] | |
|
|
3930 | 3930 | if snr0 >= snrth and snr1 >= snrth and smooth : |
|
|
3931 | 3931 | signalpn0_n0 = signalpn0 |
|
|
3932 | 3932 | signalpn0_n0[val0] = signalpn0[val0] - n0 |
|
|
3933 | mom0 = self.moments(doppler,signalpn0-n0,nProf) | |
|
|
3933 | mom0 = self.moments(doppler,signalpn0-n0,self.nProf) | |
|
|
3934 | 3934 | signalpn1_n1 = signalpn1 |
|
|
3935 | 3935 | signalpn1_n1[val1] = signalpn1[val1] - n1 |
|
|
3936 | mom1 = self.moments(doppler,signalpn1_n1,nProf) | |
|
|
3936 | mom1 = self.moments(doppler,signalpn1_n1,self.nProf) | |
|
|
3937 | 3937 | dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2. |
|
|
3938 | 3938 | dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2. |
|
|
3939 | 3939 | dataOut.data_spc = jspectra |
|
|
3940 | dataOut.spc_noise = my_noises*nProf*M | |
|
|
3941 | if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M | |
|
|
3940 | dataOut.spc_noise = my_noises*self.nProf*self.M | |
|
|
3941 | if numpy.any(proc): dataOut.spc_noise = my_noises*self.nProf*self.M | |
|
|
3942 | 3942 | if getSNR: |
|
|
3943 |
|
|
|
|
3943 | print("self.groupArray",self.groupArray.size,self.groupArray) | |
|
|
3944 | listChannels = self.groupArray.reshape((self.groupArray.size)) | |
|
|
3944 | 3945 | listChannels.sort() |
|
|
3945 | ||
|
|
3946 | dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels]) | |
|
|
3946 | # TEST | |
|
|
3947 | noise_C = numpy.zeros(self.nChannels) | |
|
|
3948 | noise_C = dataOut.getNoise() | |
|
|
3949 | print("noise_C",noise_C) | |
|
|
3950 | dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:],noise_C/(600.0*1.15))# PRUEBA *nProf*M | |
|
|
3951 | #dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise_C[listChannels])# PRUEBA *nProf*M | |
|
|
3952 | dataOut.flagNoData = False | |
|
|
3947 | 3953 | return dataOut |
|
|
3948 | 3954 | |
|
|
3949 | 3955 | def __residFunction(self, p, dp, LT, constants): |
| @@ -3956,6 +3962,11 class SpectralFitting(Operation): | |||
|
|
3956 | 3962 | |
|
|
3957 | 3963 | avg = numpy.average(z, axis=1) |
|
|
3958 | 3964 | SNR = (avg.T-noise)/noise |
|
|
3965 | print("------------------------------------------------------") | |
|
|
3966 | print("T - Noise", noise.shape) | |
|
|
3967 | print("T - Noise", 10*numpy.log10(noise[0])) | |
|
|
3968 | print("T - avg.T" , avg.T.shape) | |
|
|
3969 | print("T - avg.T" , 10*numpy.log10(avg.T[:,0])) | |
|
|
3959 | 3970 | SNR = SNR.T |
|
|
3960 | 3971 | return SNR |
|
|
3961 | 3972 | |
| @@ -4647,6 +4658,8 class EWDriftsEstimation(Operation): | |||
|
|
4647 | 4658 | |
|
|
4648 | 4659 | def run(self, dataOut, zenith, zenithCorrection,fileDrifts): |
|
|
4649 | 4660 | |
|
|
4661 | dataOut.lat=-11.95 | |
|
|
4662 | dataOut.lon=-76.87 | |
|
|
4650 | 4663 | heiRang = dataOut.heightList |
|
|
4651 | 4664 | velRadial = dataOut.data_param[:,3,:] |
|
|
4652 | 4665 | velRadialm = dataOut.data_param[:,2:4,:]*-1 |
| @@ -4721,8 +4734,9 class EWDriftsEstimation(Operation): | |||
|
|
4721 | 4734 | dataOut.heightList = heiRang1 |
|
|
4722 | 4735 | snr1 = 10*numpy.log10(SNR1[0]) |
|
|
4723 | 4736 | dataOut.data_output = winds |
|
|
4724 | #snr1 = 10*numpy.log10(SNR1[0]) | |
|
|
4737 | #snr1 = 10*numpy.log10(SNR1[0])# estaba comentado | |
|
|
4725 | 4738 | dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0])) |
|
|
4739 | print("data_snr1",dataOut.data_snr1) | |
|
|
4726 | 4740 | dataOut.utctimeInit = dataOut.utctime |
|
|
4727 | 4741 | dataOut.outputInterval = dataOut.timeInterval |
|
|
4728 | 4742 | |
General Comments 0
You need to be logged in to leave comments.
Login now