diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index a83f0f4..41ce3ea 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -26,6 +26,10 @@ from numpy import NaN from scipy.optimize.optimize import OptimizeWarning warnings.filterwarnings('ignore') +import os +import csv +from scipy import signal +import matplotlib.pyplot as plt SPEED_OF_LIGHT = 299792458 @@ -72,27 +76,47 @@ class ParametersProc(ProcessingUnit): self.dataOut.useLocalTime = self.dataIn.useLocalTime self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() + self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy() self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy() + self.dataOut.channelList = self.dataIn.channelList self.dataOut.heightList = self.dataIn.heightList + self.dataOut.ipp = self.dataIn.ipp + self.dataOut.ippSeconds = self.dataIn.ippSeconds + self.dataOut.deltaHeight = self.dataIn.deltaHeight self.dataOut.dtype = numpy.dtype([('real','= 1: + norm = normFactor[ind] + if m > 2 and m < oldfreq.size - 3: newindex = m + numpy.array([-2,-1,0,1,2]) newfreq = numpy.arange(20)/20.0*(numpy.max(freq[newindex])-numpy.min(freq[newindex]))+numpy.min(freq[newindex]) @@ -2691,6 +2738,7 @@ class SpectralMoments(Operation): power = ((spec2[valid] - n0) * fwindow[valid]).sum() fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power) + spec2 /=(norm) #compensation for sats remove snr = (spec2.mean()-n0)/n0 if (snr < 1.e-20): snr = 1.e-20 @@ -3149,6 +3197,7 @@ class SpectralFitting(Operation): self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks]) 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): + if (nicoh is None): nicoh = 1 if (graph is None): graph = 0 if (smooth is None): smooth = 0 @@ -3396,7 +3445,7 @@ class SpectralFitting(Operation): nPairs = len(crosspairs) hval=(heights >= min_hei).nonzero() ih=hval[0] - + for ih in range(hval[0][0],nHei): for ifreq in range(nProf): for ii in range(n_funct): @@ -3868,7 +3917,7 @@ class SpectralFitting(Operation): buffer3 = 0 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,Gaussian_Windowed=0): - + if not numpy.any(proc): nChannels = dataOut.nChannels nHeights= dataOut.heightList.size @@ -4058,6 +4107,8 @@ class SpectralFitting(Operation): if getSNR: listChannels = groupArray.reshape((groupArray.size)) listChannels.sort() + # norm Este factor debe estar implementado para ploteo o grabado como metadata + # norm = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels]) else: if numpy.any(taver): taver=int(taver) @@ -4307,6 +4358,8 @@ class SpectralFitting(Operation): if 0: listChannels = groupArray.reshape((groupArray.size)) listChannels.sort() + # norm Este factor debe estar implementado para ploteo o grabado como metadata + # norm = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels]) #print(dataOut.data_snr1_i) # Adding coherent echoes from possible satellites. @@ -4348,7 +4401,7 @@ class SpectralFitting(Operation): fmp=numpy.dot(LT,fm) return dp-fmp - def __getSNR(self, z, noise): + def __getSNR(self, z, noise, norm=1): avg = numpy.average(z, axis=1) SNR = (avg.T-noise)/noise @@ -4363,7 +4416,7 @@ class SpectralFitting(Operation): chisq=numpy.dot((dp-fmp).T,(dp-fmp)) return chisq -class WindProfiler_V0(Operation): +class WindProfiler(Operation): __isConfig = False @@ -4423,7 +4476,6 @@ class WindProfiler_V0(Operation): minid = listPhi.index(min(listPhi)) rango = list(range(len(phi))) - # rango = numpy.delete(rango,maxid) heiRang1 = heiRang*math.cos(phi[maxid]) heiRangAux = heiRang*math.cos(phi[minid]) @@ -4453,18 +4505,12 @@ class WindProfiler_V0(Operation): def __calculateVelUVW(self, A, velRadial): #Operacion Matricial -# velUVW = numpy.zeros((velRadial.shape[1],3)) -# for ind in range(velRadial.shape[1]): -# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind]) -# velUVW = velUVW.transpose() velUVW = numpy.zeros((A.shape[0],velRadial.shape[1])) velUVW[:,:] = numpy.dot(A,velRadial) return velUVW -# def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0): - def techniqueDBS(self, kwargs): """ Function that implements Doppler Beam Swinging (DBS) technique. @@ -4541,17 +4587,6 @@ class WindProfiler_V0(Operation): return distx, disty, dist, ang #Calculo de Matrices -# nPairs = len(pairs) -# ang1 = numpy.zeros((nPairs, 2, 1)) -# dist1 = numpy.zeros((nPairs, 2, 1)) -# -# for j in range(nPairs): -# dist1[j,0,0] = dist[pairs[j][0]] -# dist1[j,1,0] = dist[pairs[j][1]] -# ang1[j,0,0] = ang[pairs[j][0]] -# ang1[j,1,0] = ang[pairs[j][1]] -# -# return distx,disty, dist1,ang1 def __calculateVelVer(self, phase, lagTRange, _lambda): @@ -4579,26 +4614,6 @@ class WindProfiler_V0(Operation): return vel -# def __getPairsAutoCorr(self, pairsList, nChannels): -# -# pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan -# -# for l in range(len(pairsList)): -# firstChannel = pairsList[l][0] -# secondChannel = pairsList[l][1] -# -# #Obteniendo pares de Autocorrelacion -# if firstChannel == secondChannel: -# pairsAutoCorr[firstChannel] = int(l) -# -# pairsAutoCorr = pairsAutoCorr.astype(int) -# -# pairsCrossCorr = range(len(pairsList)) -# pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) -# -# return pairsAutoCorr, pairsCrossCorr - -# def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): def techniqueSA(self, kwargs): """ @@ -4626,22 +4641,10 @@ class WindProfiler_V0(Operation): _lambda = kwargs['_lambda'] #Cross Correlation pairs obtained -# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels) -# pairsArray = numpy.array(pairsList)[pairsCrossCorr] -# pairsSelArray = numpy.array(pairsSelected) -# pairs = [] -# -# #Wind estimation pairs obtained -# for i in range(pairsSelArray.shape[0]/2): -# ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0] -# ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0] -# pairs.append((ind1,ind2)) indtau = tau.shape[0]/2 tau1 = tau[:indtau,:] tau2 = tau[indtau:-1,:] -# tau1 = tau1[pairs,:] -# tau2 = tau2[pairs,:] phase1 = tau[-1,:] #--------------------------------------------------------------------- @@ -4651,11 +4654,6 @@ class WindProfiler_V0(Operation): winds = stats.nanmean(winds, axis=0) #--------------------------------------------------------------------- #Metodo General -# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth) -# #Calculo Coeficientes de Funcion de Correlacion -# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n) -# #Calculo de Velocidades -# winds = self.calculateVelUV(F,G,A,B,H) #--------------------------------------------------------------------- winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) @@ -5781,7 +5779,6 @@ class EWDriftsEstimation(Operation): rbufc=dataOut.data_paramC[:,:,0] ebufc=dataOut.data_paramC[:,:,1] - #SNR = dataOut.data_snr SNR = dataOut.data_snr1_i rbufi = dataOut.data_snr1_i velRerr = dataOut.data_error[:,4,:] @@ -7189,7 +7186,6 @@ class SMOperations(): ph0_aux = ph0 + ph1 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux)) - #First Estimation cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1)) @@ -7466,7 +7462,6 @@ class MergeProc(ProcessingUnit): #exit(1) self.dataOut.NLAG = 16 self.dataOut.NLAG = self.dataOut.data_acf.shape[1] - self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1] #print(numpy.shape(self.dataOut.data_spc)) @@ -7571,3 +7566,48 @@ class MST_Den_Conv(Operation): print("pow den shape", numpy.shape(dataOut.PowDen)) print("far den shape", numpy.shape(dataOut.FarDen)) return dataOut + +class addTxPower(Operation): + ''' + Transmited power level integrated in the dataOut ->AMISR + resolution 1 min + The power files have the pattern power_YYYYMMDD.csv + ''' + __slots__ =('isConfig','dataDatetimes','txPowers') + def __init__(self): + + Operation.__init__(self) + self.isConfig = False + self.dataDatetimes = [] + self.txPowers = [] + + def setup(self, powerFile, dutyCycle): + if not os.path.isfile(powerFile): + raise schainpy.admin.SchainError('There is no file named :{}'.format(powerFile)) + return + + with open(powerFile, newline='') as pfile: + reader = csv.reader(pfile, delimiter=',', quotechar='|') + next(reader) + for row in reader: + #'2022-10-25 00:00:00' + self.dataDatetimes.append(datetime.datetime.strptime(row[0], "%Y-%m-%d %H:%M:%S")) + self.txPowers.append(float(row[1])/dutyCycle) + self.isConfig = True + + def run(self, dataOut, path, DS=0.05): + + #dataOut.flagNoData = True + + if not(self.isConfig): + self.setup(path, DS) + + dataDate = datetime.datetime.utcfromtimestamp(dataOut.utctime).replace(second=0, microsecond=0)#no seconds + try: + indx = self.dataDatetimes.index(dataDate) + dataOut.txPower = self.txPowers[indx] + except: + log.warning("No power available for the datetime {}, setting power to 0 w", self.name) + dataOut.txPower = 0 + + return dataOut \ No newline at end of file