##// END OF EJS Templates
schain
RSS
nuevo metodo setup para atributos fijos del dataOut usado en el metodo run. Eliminacion de calculo de SNR
Alexander Valdez -
r1704:cd331037d2d0
parent child
Show More
Side by Side Unified
Context file:
r1704:cd331037d2d0 Loading diff...:
Show file before | Show file after | Hide comments
@@ -3158,10 +3158,59 class SpectralFitting(Operation):
3158 self.isConfig = False
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 self.__dataReady = False
3162 self.__dataReady = False
3163 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
3163 # new
3164 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
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 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):
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 incoh_aver[pair[1],incoh_echoes]=1
3392 incoh_aver[pair[1],incoh_echoes]=1
3344 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
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 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3395 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3348
3396
3349 nProf = dataOut.nProfiles
3397 nProf = dataOut.nProfiles
@@ -3520,6 +3568,7 class SpectralFitting(Operation):
3520 if len(valid[0]) > 0:
3568 if len(valid[0]) > 0:
3521 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3569 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3522 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3570 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3571
3523 def REM_ISOLATED_POINTS(self,array,rth):
3572 def REM_ISOLATED_POINTS(self,array,rth):
3524 if rth == None : rth = 4
3573 if rth == None : rth = 4
3525 num_prof = len(array[0,:,0])
3574 num_prof = len(array[0,:,0])
@@ -3603,37 +3652,39 class SpectralFitting(Operation):
3603 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3652 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3604 return [fmom,numpy.sqrt(smom)]
3653 return [fmom,numpy.sqrt(smom)]
3605
3654
3606
3607
3608
3609
3610 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):
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):
3656 if not self.isConfig:
3612 nChannels = dataOut.nChannels
3657 self.setup(dataOut = dataOut,groupList=groupList,path=path,file=file,filec=filec)
3613 nHeights= dataOut.heightList.size
3658 self.isConfig = True
3614 nProf = dataOut.nProfiles
3659
3615 if numpy.any(taver): taver=int(taver)
3660 if not numpy.any(proc):
3616 else : taver = 5
3661 if numpy.any(taver):
3617 tini=time.localtime(dataOut.utctime)
3662 taver = int(taver)
3663 else :
3664 taver = 5
3665 tini = time.localtime(dataOut.utctime)
3618 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3666 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3619 self.index = 0
3667 self.index = 0
3620 jspc = self.buffer
3668 jspc = self.buffer
3621 jcspc = self.buffer2
3669 jcspc = self.buffer2
3622 jnoise = self.buffer3
3670 jnoise = self.buffer3
3623 self.buffer = dataOut.data_spc
3671 self.buffer = dataOut.data_spc
3624 self.buffer2 = dataOut.data_cspc
3672 self.buffer2 = dataOut.data_cspc
3625 self.buffer3 = dataOut.noise
3673 self.buffer3 = dataOut.noise
3626 self.fint = 1
3674 self.fint = 1
3627 if numpy.any(jspc) :
3675 if numpy.any(jspc) :
3628 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
3676 jspc = numpy.reshape(jspc ,(int(len(jspc) / self.nChannels) , self.nChannels ,self.nProf,self.nHeights ))
3629 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
3677 jcspc = numpy.reshape(jcspc ,(int(len(jcspc) /int(self.nChannels/2)),int(self.nChannels/2),self.nProf,self.nHeights ))
3630 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
3678 jnoise = numpy.reshape(jnoise,(int(len(jnoise)/ self.nChannels) , self.nChannels))
3631 else:
3679 else:
3632 dataOut.flagNoData = True
3680 dataOut.flagNoData = True
3633 return dataOut
3681 return dataOut
3634 else :
3682 else :
3635 if (tini.tm_min % taver) == 0 : self.fint = 1
3683 if (tini.tm_min % taver) == 0 :
3636 else : self.fint = 0
3684 self.fint = 1
3685 else :
3686 self.fint = 0
3687
3637 self.index += 1
3688 self.index += 1
3638 if numpy.any(self.buffer):
3689 if numpy.any(self.buffer):
3639 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
3690 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
@@ -3645,85 +3696,36 class SpectralFitting(Operation):
3645 self.buffer3 = dataOut.noise
3696 self.buffer3 = dataOut.noise
3646 dataOut.flagNoData = True
3697 dataOut.flagNoData = True
3647 return dataOut
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:
3700 jnoise = jnoise/self.N# creo que falta dividirlo entre N
3699 self.isConfig = True
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 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3707 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3706 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3708 jspectra = tmp_spectra * len(jspc[:,0,0,0])
3707 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3709 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)
3710
3709 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
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)
3710 dataOut.data_spc = incoh_spectra
3712 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(self.snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3711 dataOut.data_cspc = incoh_cspectra
3713
3712 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3714 dataOut.data_spc = incoh_spectra
3713 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3715 dataOut.data_cspc = incoh_cspectra
3714 dataOut.clean_num_aver = clean_num_aver
3716 clean_num_aver = incoh_aver * len(jspc[:,0,0,0])
3715 dataOut.coh_num_aver = coh_num_aver
3717 coh_num_aver = clean_coh_aver* len(jspc[:,0,0,0])
3716 dataOut.tmp_spectra_i = incoh_spectra
3718 dataOut.clean_num_aver = clean_num_aver
3717 dataOut.tmp_cspectra_i = incoh_cspectra
3719 dataOut.coh_num_aver = coh_num_aver
3718 dataOut.tmp_spectra_c = clean_coh_spectra
3720
3719 dataOut.tmp_cspectra_c = clean_coh_cspectra
3720 #List of possible combinations
3721 #List of possible combinations
3721 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3722 #listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3722 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3723 #indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3723 if getSNR:
3724 #if getSNR:
3724 listChannels = groupArray.reshape((groupArray.size))
3725 # listChannels = groupArray.reshape((groupArray.size))
3725 listChannels.sort()
3726 # listChannels.sort()
3726 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
3727 # print("AQUI ESTOY")
3728 # dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
3727 else:
3729 else:
3728 clean_num_aver = dataOut.clean_num_aver
3730 clean_num_aver = dataOut.clean_num_aver
3729 coh_num_aver = dataOut.coh_num_aver
3731 coh_num_aver = dataOut.coh_num_aver
@@ -3762,51 +3764,49 class SpectralFitting(Operation):
3762 constants['M'] = M
3764 constants['M'] = M
3763 dataOut.constants = constants
3765 dataOut.constants = constants
3764
3766
3765 groupArray = numpy.array(groupList)
3766 dataOut.groupList = groupArray
3767 nGroups = groupArray.shape[0]
3768 #List of possible combinations
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 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3769 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3771 if dataOut.data_paramC is None:
3770 if dataOut.data_paramC is None:
3772 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
3771 dataOut.data_paramC = numpy.zeros((self.nGroups*4, self.nHeights ,2))*numpy.nan
3773 for i in range(nGroups):
3772 for i in range(self.nGroups):
3774 coord = groupArray[i,:]
3773 coord = self.groupArray[i,:]
3775 #Input data array
3774 #Input data array
3776 data = dataOut.data_spc[coord,:,:]/(M*N)
3775 data = dataOut.data_spc[coord,:,:]/(self.M*self.N)
3777 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3776 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3778
3777
3779 #Cross Spectra data array for Covariance Matrixes
3778 #Cross Spectra data array for Covariance Matrixes
3780 ind = 0
3779 ind = 0
3781 for pairs in listComb:
3780 for pairs in listComb:
3782 pairsSel = numpy.array([coord[x],coord[y]])
3781 pairsSel = numpy.array([coord[x],coord[y]])
3783 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
3782 indCross[ind] = int(numpy.where(numpy.all(self.pairsArray == pairsSel, axis = 1))[0][0])
3784 ind += 1
3783 ind += 1
3785 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
3784 dataCross = dataOut.data_cspc[indCross,:,:]/(self.M*self.N)
3786 dataCross = dataCross**2
3785 dataCross = dataCross**2
3787 nhei = nHeights
3786 nhei = self.nHeights
3788 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
3787 poweri = numpy.sum(dataOut.data_spc[:,1:self.nProf-0,:],axis=1)/clean_num_aver[:,:]
3788
3789 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3789 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3790 n0i = numpy.nanmin(poweri[0+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])/(nProf-1)
3791 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(self.nProf-1)
3792 n0 = n0i
3792 n0 = n0i
3793 n1= n1i
3793 n1= n1i
3794 my_noises[2*i+0] = n0
3794 my_noises[2*i+0] = n0
3795 my_noises[2*i+1] = n1
3795 my_noises[2*i+1] = n1
3796 snrth = -15.0 # -4 -16 -25
3796 snrth = -15.0 # -4 -16 -25
3797 snrth = 10**(snrth/10.0)
3797 snrth = 10**(snrth/10.0)
3798 jvelr = numpy.zeros(nHeights, dtype = 'float')
3798 jvelr = numpy.zeros(self.nHeights, dtype = 'float')
3799 hvalid = [0]
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 smooth = clean_num_aver[i+1,h]
3803 smooth = clean_num_aver[i+1,h]
3804 signalpn0 = (dataOut.data_spc[i*2,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:(nProf-0),h])/smooth
3805 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3806 signal0 = signalpn0-n0
3806 signal0 = signalpn0-n0
3807 signal1 = signalpn1-n1
3807 signal1 = signalpn1-n1
3808 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3808 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3809 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3809 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3810 gamma = coh2[:,h]
3810 gamma = coh2[:,h]
3811 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3811 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3812 if len(indxs) >0:
3812 if len(indxs) >0:
@@ -3817,22 +3817,22 class SpectralFitting(Operation):
3817 else:
3817 else:
3818 maxp0 = numpy.argmax(signal0)
3818 maxp0 = numpy.argmax(signal0)
3819 maxp1 = numpy.argmax(signal1)
3819 maxp1 = numpy.argmax(signal1)
3820 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
3820 jvelr[h] = (self.absc[maxp0]+self.absc[maxp1])/2.
3821 else: jvelr[h] = absc[0]
3821 else: jvelr[h] = self.absc[0]
3822 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3822 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3823 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3823 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3824
3824
3825 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3825 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3826 else: fd0 = numpy.nan
3826 else: fd0 = numpy.nan
3827 for h in range(nHeights):
3827 for h in range(self.nHeights):
3828 d = data[:,h]
3828 d = data[:,h]
3829 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
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
3830 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3831 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3831 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3832 signal0 = signalpn0-n0
3832 signal0 = signalpn0-n0
3833 signal1 = signalpn1-n1
3833 signal1 = signalpn1-n1
3834 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3834 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3835 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3835 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3836 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3836 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3837 #Covariance Matrix
3837 #Covariance Matrix
3838 D = numpy.diag(d**2)
3838 D = numpy.diag(d**2)
@@ -3867,7 +3867,7 class SpectralFitting(Operation):
3867 if (h>6)and(error1[3]<25):
3867 if (h>6)and(error1[3]<25):
3868 p0 = dataOut.data_param[i,:,h-1].copy()
3868 p0 = dataOut.data_param[i,:,h-1].copy()
3869 else:
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 p0[3] = fd0
3871 p0[3] = fd0
3872
3872
3873 if filec != None:
3873 if filec != None:
@@ -3875,12 +3875,12 class SpectralFitting(Operation):
3875
3875
3876 try:
3876 try:
3877 #Least Squares
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 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3879 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3880 #Chi square error
3880 #Chi square error
3881 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
3881 error0 = numpy.sum(infodict['fvec']**2)/(2*self.N)
3882 #Error with Jacobian
3882 #Error with Jacobian
3883 error1 = self.library.errorFunction(minp,constants,LT)
3883 error1 = self.library.errorFunction(minp,self.constants,LT)
3884
3884
3885 except:
3885 except:
3886 minp = p0*numpy.nan
3886 minp = p0*numpy.nan
@@ -3888,29 +3888,29 class SpectralFitting(Operation):
3888 error1 = p0*numpy.nan
3888 error1 = p0*numpy.nan
3889 else :
3889 else :
3890 data_spc = dataOut.data_spc[coord,:,h]
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 minp = p0*numpy.nan
3892 minp = p0*numpy.nan
3893 error0 = numpy.nan
3893 error0 = numpy.nan
3894 error1 = p0*numpy.nan
3894 error1 = p0*numpy.nan
3895 if dataOut.data_param is None:
3895 if dataOut.data_param is None:
3896 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
3896 dataOut.data_param = numpy.zeros((self.nGroups, p0.size, self.nHeights ))*numpy.nan
3897 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
3897 dataOut.data_error = numpy.zeros((self.nGroups, p0.size + 1, self.nHeights ))*numpy.nan
3898 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3898 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3899 dataOut.data_param[i,:,h] = minp
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 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3902 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3903 dataOut.data_paramC[4*i,ht,1] = smooth
3903 dataOut.data_paramC[4*i,ht,1] = smooth
3904 signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
3904 signalpn0 = (clean_coh_spectra[i*2 ,1:(self.nProf-0),ht])/smooth #coh_spectra
3905 signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
3905 signalpn1 = (clean_coh_spectra[i*2+1,1:(self.nProf-0),ht])/smooth
3906 val0 = (signalpn0 > 0).nonzero()
3906 val0 = (signalpn0 > 0).nonzero()
3907 val0 = val0[0]
3907 val0 = val0[0]
3908 if len(val0) == 0 : val0_npoints = nProf
3908 if len(val0) == 0 : val0_npoints = self.nProf
3909 else : val0_npoints = len(val0)
3909 else : val0_npoints = len(val0)
3910
3910
3911 val1 = (signalpn1 > 0).nonzero()
3911 val1 = (signalpn1 > 0).nonzero()
3912 val1 = val1[0]
3912 val1 = val1[0]
3913 if len(val1) == 0 : val1_npoints = nProf
3913 if len(val1) == 0 : val1_npoints = self.nProf
3914 else : val1_npoints = len(val1)
3914 else : val1_npoints = len(val1)
3915
3915
3916 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
3916 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
@@ -3924,26 +3924,32 class SpectralFitting(Operation):
3924 vali = (signal1 < 0).nonzero()
3924 vali = (signal1 < 0).nonzero()
3925 vali = vali[0]
3925 vali = vali[0]
3926 if len(vali) > 0 : signal1[vali] = 0
3926 if len(vali) > 0 : signal1[vali] = 0
3927 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3927 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3928 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3928 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3929 doppler = absc[1:]
3929 doppler = self.absc[1:]
3930 if snr0 >= snrth and snr1 >= snrth and smooth :
3930 if snr0 >= snrth and snr1 >= snrth and smooth :
3931 signalpn0_n0 = signalpn0
3931 signalpn0_n0 = signalpn0
3932 signalpn0_n0[val0] = signalpn0[val0] - n0
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 signalpn1_n1 = signalpn1
3934 signalpn1_n1 = signalpn1
3935 signalpn1_n1[val1] = signalpn1[val1] - n1
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 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3937 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3938 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
3938 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
3939 dataOut.data_spc = jspectra
3939 dataOut.data_spc = jspectra
3940 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*nProf*M
3941 if numpy.any(proc): dataOut.spc_noise = my_noises*self.nProf*self.M
3942 if getSNR:
3942 if getSNR:
3943 listChannels = groupArray.reshape((groupArray.size))
3943 print("self.groupArray",self.groupArray.size,self.groupArray)
3944 listChannels = self.groupArray.reshape((self.groupArray.size))
3944 listChannels.sort()
3945 listChannels.sort()
3945
3946 # TEST
3946 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
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 return dataOut
3953 return dataOut
3948
3954
3949 def __residFunction(self, p, dp, LT, constants):
3955 def __residFunction(self, p, dp, LT, constants):
@@ -3956,6 +3962,11 class SpectralFitting(Operation):
3956
3962
3957 avg = numpy.average(z, axis=1)
3963 avg = numpy.average(z, axis=1)
3958 SNR = (avg.T-noise)/noise
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 SNR = SNR.T
3970 SNR = SNR.T
3960 return SNR
3971 return SNR
3961
3972
@@ -4647,6 +4658,8 class EWDriftsEstimation(Operation):
4647
4658
4648 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4659 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4649
4660
4661 dataOut.lat=-11.95
4662 dataOut.lon=-76.87
4650 heiRang = dataOut.heightList
4663 heiRang = dataOut.heightList
4651 velRadial = dataOut.data_param[:,3,:]
4664 velRadial = dataOut.data_param[:,3,:]
4652 velRadialm = dataOut.data_param[:,2:4,:]*-1
4665 velRadialm = dataOut.data_param[:,2:4,:]*-1
@@ -4721,8 +4734,9 class EWDriftsEstimation(Operation):
4721 dataOut.heightList = heiRang1
4734 dataOut.heightList = heiRang1
4722 snr1 = 10*numpy.log10(SNR1[0])
4735 snr1 = 10*numpy.log10(SNR1[0])
4723 dataOut.data_output = winds
4736 dataOut.data_output = winds
4724 #snr1 = 10*numpy.log10(SNR1[0])
4737 #snr1 = 10*numpy.log10(SNR1[0])# estaba comentado
4725 dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
4738 dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
4739 print("data_snr1",dataOut.data_snr1)
4726 dataOut.utctimeInit = dataOut.utctime
4740 dataOut.utctimeInit = dataOut.utctime
4727 dataOut.outputInterval = dataOut.timeInterval
4741 dataOut.outputInterval = dataOut.timeInterval
4728
4742
General Comments 0
You need to be logged in to leave comments. Login now