From ae31287a1940c7fddbee1f597225a950278e970d 2016-09-19 15:51:33 From: Julio Valdez Date: 2016-09-19 15:51:33 Subject: [PATCH] Parameters libraries reorganized, now each operation is a separated class --- diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 1213128..2278caa 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -111,39 +111,48 @@ class ParametersProc(ProcessingUnit): self.__updateObjFromInput() self.dataOut.utctimeInit = self.dataIn.utctime self.dataOut.paramInterval = self.dataIn.timeInterval + + return - #------------------- Get Moments ---------------------------------- - def GetMoments(self, channelList = None): - ''' - Function GetMoments() +class SpectralMoments(Operation): + + ''' + Function SpectralMoments() + + Calculates moments (power, mean, standard deviation) and SNR of the signal + + Type of dataIn: Spectra Input: channelList : simple channel list to select e.g. [2,3,7] - self.dataOut.data_pre - self.dataOut.abscissaList - self.dataOut.noise + self.dataOut.data_pre : Spectral data + self.dataOut.abscissaList : List of frequencies + self.dataOut.noise : Noise level per channel Affected: - self.dataOut.data_param - self.dataOut.data_SNR + self.dataOut.data_param : Parameters per channel + self.dataOut.data_SNR : SNR per channel - ''' - self.dataOut.data_pre = self.dataOut.data_pre[0] - data = self.dataOut.data_pre - absc = self.dataOut.abscissaList[:-1] - noise = self.dataOut.noise + ''' + + def run(self, dataOut, channelList=None): + + dataOut.data_pre = dataOut.data_pre[0] + data = dataOut.data_pre + absc = dataOut.abscissaList[:-1] + noise = dataOut.noise data_param = numpy.zeros((data.shape[0], 4, data.shape[2])) if channelList== None: - channelList = self.dataIn.channelList - self.dataOut.channelList = channelList + channelList = dataOut.channelList + dataOut.channelList = channelList for ind in channelList: data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind]) - self.dataOut.data_param = data_param[:,1:,:] - self.dataOut.data_SNR = data_param[:,0] + dataOut.data_param = data_param[:,1:,:] + dataOut.data_SNR = data_param[:,0] return 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): @@ -221,6 +230,7 @@ class ParametersProc(ProcessingUnit): #------------------ Get SA Parameters -------------------------- def GetSAParameters(self): + #SA en frecuencia pairslist = self.dataOut.groupList num_pairs = len(pairslist) @@ -249,61 +259,64 @@ class ParametersProc(ProcessingUnit): cspc_par[i,:,:] = __calculateMoments(norm_cspectra) #------------------- Get Lags ---------------------------------- - def GetLags(self): - ''' - Function GetMoments() +class SALags(Operation): + ''' + Function GetMoments() - Input: - self.dataOut.data_pre - self.dataOut.abscissaList - self.dataOut.noise - self.dataOut.normFactor - self.dataOut.data_SNR - self.dataOut.groupList - self.dataOut.nChannels - - Affected: - self.dataOut.data_param - - ''' - - data = self.dataOut.data_pre - normFactor = self.dataOut.normFactor - nHeights = self.dataOut.nHeights - absc = self.dataOut.abscissaList[:-1] - noise = self.dataOut.noise - SNR = self.dataOut.data_SNR - pairsList = self.dataOut.groupList - nChannels = self.dataOut.nChannels - pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) - self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights)) + Input: + self.dataOut.data_pre + self.dataOut.abscissaList + self.dataOut.noise + self.dataOut.normFactor + self.dataOut.data_SNR + self.dataOut.groupList + self.dataOut.nChannels + + Affected: + self.dataOut.data_param + + ''' + def run(self, dataOut): + data = dataOut.data_pre + data_acf = dataOut.data_pre[0] + data_ccf = dataOut.data_pre[1] + + normFactor = dataOut.normFactor + nHeights = dataOut.nHeights + absc = dataOut.abscissaList[:-1] + noise = dataOut.noise + SNR = dataOut.data_SNR +# pairsList = dataOut.groupList + nChannels = dataOut.nChannels +# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) + dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights)) dataNorm = numpy.abs(data) for l in range(len(pairsList)): dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:] self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc) - self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc) + self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc) return - 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 __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 __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange): @@ -331,1016 +344,913 @@ class ParametersProc(ProcessingUnit): return tau - def __calculateLag1Phase(self, data, pairs, lagTRange): - data1 = stats.nanmean(data[pairs,:,:], axis = 0) + def __calculateLag1Phase(self, data, lagTRange): + data1 = stats.nanmean(data, axis = 0) lag1 = numpy.where(lagTRange == 0)[0][0] + 1 phase = numpy.angle(data1[lag1,:]) return phase - #------------------- Detect Meteors ------------------------------ - def MeteorDetection(self, hei_ref = None, tauindex = 0, - phaseOffsets = None, - cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, - noise_timeStep = 4, noise_multiple = 4, - multDet_timeLimit = 1, multDet_rangeLimit = 3, - phaseThresh = 20, SNRThresh = 5, - hmin = 50, hmax=150, azimuth = 0, -# channel25X = 0, channel20X = 4, channelCentX = 2, -# channel25Y = 1, channel20Y = 3, channelCentY = 2, - channelPositions = None) : +class SpectralFitting(Operation): + ''' + Function GetMoments() - ''' - Function DetectMeteors() - Project developed with paper: - HOLDSWORTH ET AL. 2004 - Input: - self.dataOut.data_pre - - centerReceiverIndex: From the channels, which is the center receiver - - hei_ref: Height reference for the Beacon signal extraction - tauindex: - predefinedPhaseShifts: Predefined phase offset for the voltge signals - - cohDetection: Whether to user Coherent detection or not - cohDet_timeStep: Coherent Detection calculation time step - cohDet_thresh: Coherent Detection phase threshold to correct phases - - noise_timeStep: Noise calculation time step - noise_multiple: Noise multiple to define signal threshold - - multDet_timeLimit: Multiple Detection Removal time limit in seconds - multDet_rangeLimit: Multiple Detection Removal range limit in km - - phaseThresh: Maximum phase difference between receiver to be consider a meteor - SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor - - hmin: Minimum Height of the meteor to use it in the further wind estimations - hmax: Maximum Height of the meteor to use it in the further wind estimations - azimuth: Azimuth angle correction - - Affected: - self.dataOut.data_param - - Rejection Criteria (Errors): - 0: No error; analysis OK - 1: SNR < SNR threshold - 2: angle of arrival (AOA) ambiguously determined - 3: AOA estimate not feasible - 4: Large difference in AOAs obtained from different antenna baselines - 5: echo at start or end of time series - 6: echo less than 5 examples long; too short for analysis - 7: echo rise exceeds 0.3s - 8: echo decay time less than twice rise time - 9: large power level before echo - 10: large power level after echo - 11: poor fit to amplitude for estimation of decay time - 12: poor fit to CCF phase variation for estimation of radial drift velocity - 13: height unresolvable echo: not valid height within 70 to 110 km - 14: height ambiguous echo: more then one possible height within 70 to 110 km - 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s - 16: oscilatory echo, indicating event most likely not an underdense echo - - 17: phase difference in meteor Reestimation - - Data Storage: - Meteors for Wind Estimation (8): - Utc Time | Range Height - Azimuth Zenith errorCosDir - VelRad errorVelRad - Phase0 Phase1 Phase2 Phase3 - TypeError - - ''' - #Getting Pairslist - if channelPositions == None: -# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T - channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella - meteorOps = MeteorOperations() - pairslist0, distances = meteorOps.getPhasePairs(channelPositions) + Output: + Variables modified: + ''' + + def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None): - #Get Beacon signal - newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex]) - if hei_ref != None: - newheis = numpy.where(self.dataOut.heightList>hei_ref) + if path != None: + sys.path.append(path) + self.dataOut.library = importlib.import_module(file) - heiRang = self.dataOut.getHeiRange() - -# nChannel = self.dataOut.nChannels -# for i in range(nChannel): -# if i != centerReceiverIndex: -# pairslist.append((centerReceiverIndex,i)) - - #****************REMOVING HARDWARE PHASE DIFFERENCES*************** - # see if the user put in pre defined phase shifts - voltsPShift = self.dataOut.data_pre.copy() + #To be inserted as a parameter + groupArray = numpy.array(groupList) +# groupArray = numpy.array([[0,1],[2,3]]) + self.dataOut.groupList = groupArray -# if predefinedPhaseShifts != None: -# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 -# -# # elif beaconPhaseShifts: -# # #get hardware phase shifts using beacon signal -# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) -# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) -# -# else: -# hardwarePhaseShifts = numpy.zeros(5) -# -# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex') -# for i in range(self.dataOut.data_pre.shape[0]): -# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i]) - - #******************END OF REMOVING HARDWARE PHASE DIFFERENCES********* - - #Remove DC - voltsDC = numpy.mean(voltsPShift,1) - voltsDC = numpy.mean(voltsDC,1) - for i in range(voltsDC.shape[0]): - voltsPShift[i] = voltsPShift[i] - voltsDC[i] - - #Don't considerate last heights, theyre used to calculate Hardware Phase Shift - voltsPShift = voltsPShift[:,:,:newheis[0][0]] + nGroups = groupArray.shape[0] + nChannels = self.dataIn.nChannels + nHeights=self.dataIn.heightList.size - #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) ********** - #Coherent Detection - if cohDetection: - #use coherent detection to get the net power - cohDet_thresh = cohDet_thresh*numpy.pi/180 - voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh) + #Parameters Array + self.dataOut.data_param = None - #Non-coherent detection! - powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0) - #********** END OF COH/NON-COH POWER CALCULATION********************** - - #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS **************** - #Get noise - noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval) -# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval) - #Get signal threshold - signalThresh = noise_multiple*noise - #Meteor echoes detection - listMeteors = self.__findMeteors(powerNet, signalThresh) - #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********** + #Set constants + constants = self.dataOut.library.setConstants(self.dataIn) + self.dataOut.constants = constants + M = self.dataIn.normFactor + N = self.dataIn.nFFTPoints + ippSeconds = self.dataIn.ippSeconds + K = self.dataIn.nIncohInt + pairsArray = numpy.array(self.dataIn.pairsList) - #************** REMOVE MULTIPLE DETECTIONS (3.5) *************************** - #Parameters - heiRange = self.dataOut.getHeiRange() - rangeInterval = heiRange[1] - heiRange[0] - rangeLimit = multDet_rangeLimit/rangeInterval - timeLimit = multDet_timeLimit/self.dataOut.timeInterval - #Multiple detection removals - listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit) - #************ END OF REMOVE MULTIPLE DETECTIONS ********************** + #List of possible combinations + listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) + indCross = numpy.zeros(len(list(listComb)), dtype = 'int') - #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ******************** - #Parameters - phaseThresh = phaseThresh*numpy.pi/180 - thresh = [phaseThresh, noise_multiple, SNRThresh] - #Meteor reestimation (Errors N 1, 6, 12, 17) - listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency) -# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise) - #Estimation of decay times (Errors N 7, 8, 11) - listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency) - #******************* END OF METEOR REESTIMATION ******************* + if getSNR: + listChannels = groupArray.reshape((groupArray.size)) + listChannels.sort() + noise = self.dataIn.getNoise() + self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) - #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) ************************** - #Calculating Radial Velocity (Error N 15) - radialStdThresh = 10 - listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval) - - if len(listMeteors4) > 0: - #Setting New Array - date = self.dataOut.utctime - arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) - - #Correcting phase offset - if phaseOffsets != None: - phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 - arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) + for i in range(nGroups): + coord = groupArray[i,:] - #Second Pairslist - pairsList = [] - pairx = (0,1) - pairy = (2,3) - pairsList.append(pairx) - pairsList.append(pairy) + #Input data array + data = self.dataIn.data_spc[coord,:,:]/(M*N) + data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) - jph = numpy.array([0,0,0,0]) - h = (hmin,hmax) - arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) + #Cross Spectra data array for Covariance Matrixes + ind = 0 + for pairs in listComb: + pairsSel = numpy.array([coord[x],coord[y]]) + indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) + ind += 1 + dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N) + dataCross = dataCross**2/K -# #Calculate AOA (Error N 3, 4) -# #JONES ET AL. 1998 -# error = arrayParameters[:,-1] -# AOAthresh = numpy.pi/8 -# phases = -arrayParameters[:,9:13] -# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) -# -# #Calculate Heights (Error N 13 and 14) -# error = arrayParameters[:,-1] -# Ranges = arrayParameters[:,2] -# zenith = arrayParameters[:,5] -# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) -# error = arrayParameters[:,-1] - #********************* END OF PARAMETERS CALCULATION ************************** + for h in range(nHeights): +# print self.dataOut.heightList[h] - #***************************+ PASS DATA TO NEXT STEP ********************** -# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) - self.dataOut.data_param = arrayParameters - - if arrayParameters == None: - self.dataOut.flagNoData = True - else: - self.dataOut.flagNoData = True - + #Input + d = data[:,h] + + #Covariance Matrix + D = numpy.diag(d**2/K) + ind = 0 + for pairs in listComb: + #Coordinates in Covariance Matrix + x = pairs[0] + y = pairs[1] + #Channel Index + S12 = dataCross[ind,:,h] + D12 = numpy.diag(S12) + #Completing Covariance Matrix with Cross Spectras + D[x*N:(x+1)*N,y*N:(y+1)*N] = D12 + D[y*N:(y+1)*N,x*N:(x+1)*N] = D12 + ind += 1 + Dinv=numpy.linalg.inv(D) + L=numpy.linalg.cholesky(Dinv) + LT=L.T + + dp = numpy.dot(LT,d) + + #Initial values + data_spc = self.dataIn.data_spc[coord,:,h] + + if (h>0)and(error1[3]<5): + p0 = self.dataOut.data_param[i,:,h-1] + else: + p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) + + try: + #Least Squares + minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) +# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) + #Chi square error + error0 = numpy.sum(infodict['fvec']**2)/(2*N) + #Error with Jacobian + error1 = self.dataOut.library.errorFunction(minp,constants,LT) + except: + minp = p0*numpy.nan + error0 = numpy.nan + error1 = p0*numpy.nan + + #Save + if self.dataOut.data_param == None: + self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan + self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan + + self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) + self.dataOut.data_param[i,:,h] = minp return - def CorrectMeteorPhases(self, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None): - - arrayParameters = self.dataOut.data_param - pairsList = [] - pairx = (0,1) - pairy = (2,3) - pairsList.append(pairx) - pairsList.append(pairy) - jph = numpy.zeros(4) - - phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 -# arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) - arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets))) - - meteorOps = MeteorOperations() - if channelPositions == None: -# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T - channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella - - pairslist0, distances = meteorOps.getPhasePairs(channelPositions) - h = (hmin,hmax) + def __residFunction(self, p, dp, LT, constants): - arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) + fm = self.dataOut.library.modelFunction(p, constants) + fmp=numpy.dot(LT,fm) + + return dp-fmp + + def __getSNR(self, z, noise): - self.dataOut.data_param = arrayParameters - return + avg = numpy.average(z, axis=1) + SNR = (avg.T-noise)/noise + SNR = SNR.T + return SNR - def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n): - - minIndex = min(newheis[0]) - maxIndex = max(newheis[0]) - - voltage = voltage0[:,:,minIndex:maxIndex+1] - nLength = voltage.shape[1]/n - nMin = 0 - nMax = 0 - phaseOffset = numpy.zeros((len(pairslist),n)) - - for i in range(n): - nMax += nLength - phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0])) - phaseCCF = numpy.mean(phaseCCF, axis = 2) - phaseOffset[:,i] = phaseCCF.transpose() - nMin = nMax -# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist) - - #Remove Outliers - factor = 2 - wt = phaseOffset - signal.medfilt(phaseOffset,(1,5)) - dw = numpy.std(wt,axis = 1) - dw = dw.reshape((dw.size,1)) - ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) - phaseOffset[ind] = numpy.nan - phaseOffset = stats.nanmean(phaseOffset, axis=1) + def __chisq(p,chindex,hindex): + #similar to Resid but calculates CHI**2 + [LT,d,fm]=setupLTdfm(p,chindex,hindex) + dp=numpy.dot(LT,d) + fmp=numpy.dot(LT,fm) + chisq=numpy.dot((dp-fmp).T,(dp-fmp)) + return chisq + +class WindProfiler(Operation): + + __isConfig = False - return phaseOffset + __initime = None + __lastdatatime = None + __integrationtime = None - def __shiftPhase(self, data, phaseShift): - #this will shift the phase of a complex number - dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j) - return dataShifted + __buffer = None - def __estimatePhaseDifference(self, array, pairslist): - nChannel = array.shape[0] - nHeights = array.shape[2] - numPairs = len(pairslist) -# phaseCCF = numpy.zeros((nChannel, 5, nHeights)) - phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2])) + __dataReady = False + + __firstdata = None + + n = None + + def __init__(self): + Operation.__init__(self) + + def __calculateCosDir(self, elev, azim): + zen = (90 - elev)*numpy.pi/180 + azim = azim*numpy.pi/180 + cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) + cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) - #Correct phases - derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:] - indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) + signX = numpy.sign(numpy.cos(azim)) + signY = numpy.sign(numpy.sin(azim)) - if indDer[0].shape[0] > 0: - for i in range(indDer[0].shape[0]): - signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]]) - phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi + cosDirX = numpy.copysign(cosDirX, signX) + cosDirY = numpy.copysign(cosDirY, signY) + return cosDirX, cosDirY -# for j in range(numSides): -# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2]) -# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux) -# - #Linear - phaseInt = numpy.zeros((numPairs,1)) - angAllCCF = phaseCCF[:,[0,1,3,4],0] - for j in range(numPairs): - fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:]) - phaseInt[j] = fit[1] - #Phase Differences - phaseDiff = phaseInt - phaseCCF[:,2,:] - phaseArrival = phaseInt.reshape(phaseInt.size) + def __calculateAngles(self, theta_x, theta_y, azimuth): + + dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) + zenith_arr = numpy.arccos(dir_cosw) + azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 - #Dealias - phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival)) -# indAlias = numpy.where(phaseArrival > numpy.pi) -# phaseArrival[indAlias] -= 2*numpy.pi -# indAlias = numpy.where(phaseArrival < -numpy.pi) -# phaseArrival[indAlias] += 2*numpy.pi + dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) + dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) - return phaseDiff, phaseArrival - - def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh): - #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power - #find the phase shifts of each channel over 1 second intervals - #only look at ranges below the beacon signal - numProfPerBlock = numpy.ceil(timeSegment/timeInterval) - numBlocks = int(volts.shape[1]/numProfPerBlock) - numHeights = volts.shape[2] - nChannel = volts.shape[0] - voltsCohDet = volts.copy() + return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw + + def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): - pairsarray = numpy.array(pairslist) - indSides = pairsarray[:,1] -# indSides = numpy.array(range(nChannel)) -# indSides = numpy.delete(indSides, indCenter) # -# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0) - listBlocks = numpy.array_split(volts, numBlocks, 1) + if horOnly: + A = numpy.c_[dir_cosu,dir_cosv] + else: + A = numpy.c_[dir_cosu,dir_cosv,dir_cosw] + A = numpy.asmatrix(A) + A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose() + + return A1 + + def __correctValues(self, heiRang, phi, velRadial, SNR): + listPhi = phi.tolist() + maxid = listPhi.index(max(listPhi)) + minid = listPhi.index(min(listPhi)) - startInd = 0 - endInd = 0 + rango = range(len(phi)) + # rango = numpy.delete(rango,maxid) + + heiRang1 = heiRang*math.cos(phi[maxid]) + heiRangAux = heiRang*math.cos(phi[minid]) + indOut = (heiRang1 < heiRangAux[0]).nonzero() + heiRang1 = numpy.delete(heiRang1,indOut) + + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) + SNR1 = numpy.zeros([len(phi),len(heiRang1)]) + + for i in rango: + x = heiRang*math.cos(phi[i]) + y1 = velRadial[i,:] + f1 = interpolate.interp1d(x,y1,kind = 'cubic') - for i in range(numBlocks): - startInd = endInd - endInd = endInd + listBlocks[i].shape[1] + x1 = heiRang1 + y11 = f1(x1) - arrayBlock = listBlocks[i] -# arrayBlockCenter = listCenter[i] + y2 = SNR[i,:] + f2 = interpolate.interp1d(x,y2,kind = 'cubic') + y21 = f2(x1) - #Estimate the Phase Difference - phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist) - #Phase Difference RMS - arrayPhaseRMS = numpy.abs(phaseDiff) - phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0) - indPhase = numpy.where(phaseRMSaux==4) - #Shifting - if indPhase[0].shape[0] > 0: - for j in range(indSides.size): - arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose()) - voltsCohDet[:,startInd:endInd,:] = arrayBlock - - return voltsCohDet - - def __calculateCCF(self, volts, pairslist ,laglist): + velRadial1[i,:] = y11 + SNR1[i,:] = y21 + + return heiRang1, velRadial1, SNR1 + + def __calculateVelUVW(self, A, velRadial): - nHeights = volts.shape[2] - nPoints = volts.shape[1] - voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex') + #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) - for i in range(len(pairslist)): - volts1 = volts[pairslist[i][0]] - volts2 = volts[pairslist[i][1]] - - for t in range(len(laglist)): - idxT = laglist[t] - if idxT >= 0: - vStacked = numpy.vstack((volts2[idxT:,:], - numpy.zeros((idxT, nHeights),dtype='complex'))) - else: - vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'), - volts2[:(nPoints + idxT),:])) - voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0) + + return velUVW - vStacked = None - return voltsCCF +# def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0): + def techniqueDBS(self, velRadial0, heiRang, SNR0, kwargs): + """ + Function that implements Doppler Beam Swinging (DBS) technique. - def __getNoise(self, power, timeSegment, timeInterval): - numProfPerBlock = numpy.ceil(timeSegment/timeInterval) - numBlocks = int(power.shape[0]/numProfPerBlock) - numHeights = power.shape[1] - - listPower = numpy.array_split(power, numBlocks, 0) - noise = numpy.zeros((power.shape[0], power.shape[1])) - noise1 = numpy.zeros((power.shape[0], power.shape[1])) + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR - startInd = 0 - endInd = 0 + Output: Winds estimation (Zonal, Meridional and Vertical) - for i in range(numBlocks): #split por canal - startInd = endInd - endInd = endInd + listPower[i].shape[0] - - arrayBlock = listPower[i] - noiseAux = numpy.mean(arrayBlock, 0) -# noiseAux = numpy.median(noiseAux) -# noiseAux = numpy.mean(arrayBlock) - noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux - - noiseAux1 = numpy.mean(arrayBlock) - noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1 - - return noise, noise1 - - def __findMeteors(self, power, thresh): - nProf = power.shape[0] - nHeights = power.shape[1] - listMeteors = [] + Parameters affected: Winds, height range, SNR + """ - for i in range(nHeights): - powerAux = power[:,i] - threshAux = thresh[:,i] - - indUPthresh = numpy.where(powerAux > threshAux)[0] - indDNthresh = numpy.where(powerAux <= threshAux)[0] - - j = 0 - - while (j < indUPthresh.size - 2): - if (indUPthresh[j + 2] == indUPthresh[j] + 2): - indDNAux = numpy.where(indDNthresh > indUPthresh[j]) - indDNthresh = indDNthresh[indDNAux] - - if (indDNthresh.size > 0): - indEnd = indDNthresh[0] - 1 - indInit = indUPthresh[j] - - meteor = powerAux[indInit:indEnd + 1] - indPeak = meteor.argmax() + indInit - FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0))) - - listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!! - j = numpy.where(indUPthresh == indEnd)[0] + 1 - else: j+=1 - else: j+=1 - - return listMeteors + if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'): + theta_x = numpy.array(kwargs['dirCosx']) + theta_y = numpy.array(kwargs['dirCosy']) + else: + elev = numpy.array(kwargs['elevation']) + azim = numpy.array(kwargs['azimuth']) + theta_x, theta_y = self.__calculateCosDir(elev, azim) + azimuth = kwargs['correctAzimuth'] + if kwargs.has_key('horizontalOnly'): + horizontalOnly = kwargs['horizontalOnly'] + else: horizontalOnly = False + if kwargs.has_key('correctFactor'): + correctFactor = kwargs['correctFactor'] + else: correctFactor = 1 + if kwargs.has_key('channelList'): + channelList = kwargs['channelList'] + if len(channelList) == 2: + horizontalOnly = True + arrayChannel = numpy.array(channelList) + param = param[arrayChannel,:,:] + theta_x = theta_x[arrayChannel] + theta_y = theta_y[arrayChannel] + + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth) + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0) + A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) + + #Calculo de Componentes de la velocidad con DBS + winds = self.__calculateVelUVW(A,velRadial1) + + return winds, heiRang1, SNR1 - def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit): + def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None): - arrayMeteors = numpy.asarray(listMeteors) - listMeteors1 = [] + posx = numpy.asarray(posx) + posy = numpy.asarray(posy) - while arrayMeteors.shape[0] > 0: - FLAs = arrayMeteors[:,4] - maxFLA = FLAs.argmax() - listMeteors1.append(arrayMeteors[maxFLA,:]) - - MeteorInitTime = arrayMeteors[maxFLA,1] - MeteorEndTime = arrayMeteors[maxFLA,3] - MeteorHeight = arrayMeteors[maxFLA,0] - - #Check neighborhood - maxHeightIndex = MeteorHeight + rangeLimit - minHeightIndex = MeteorHeight - rangeLimit - minTimeIndex = MeteorInitTime - timeLimit - maxTimeIndex = MeteorEndTime + timeLimit - - #Check Heights - indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex) - indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex) - indBoth = numpy.where(numpy.logical_and(indTime,indHeight)) + #Rotacion Inversa para alinear con el azimuth + if azimuth!= None: + azimuth = azimuth*math.pi/180 + posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth) + posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth) + else: + posx1 = posx + posy1 = posy + + #Calculo de Distancias + distx = numpy.zeros(pairsCrossCorr.size) + disty = numpy.zeros(pairsCrossCorr.size) + dist = numpy.zeros(pairsCrossCorr.size) + ang = numpy.zeros(pairsCrossCorr.size) + + for i in range(pairsCrossCorr.size): + distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]] + disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]] + dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) + ang[i] = numpy.arctan2(disty[i],distx[i]) + #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]] - arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0) + return distx,disty, dist1,ang1 + + def __calculateVelVer(self, phase, lagTRange, _lambda): + + Ts = lagTRange[1] - lagTRange[0] + velW = -_lambda*phase/(4*math.pi*Ts) - return listMeteors1 + return velW - def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency): - numHeights = volts.shape[2] - nChannel = volts.shape[0] + def __calculateVelHorDir(self, dist, tau1, tau2, ang): + nPairs = tau1.shape[0] + vel = numpy.zeros((nPairs,3,tau1.shape[2])) - thresholdPhase = thresh[0] - thresholdNoise = thresh[1] - thresholdDB = float(thresh[2]) + angCos = numpy.cos(ang) + angSin = numpy.sin(ang) - thresholdDB1 = 10**(thresholdDB/10) - pairsarray = numpy.array(pairslist) - indSides = pairsarray[:,1] + vel0 = dist*tau1/(2*tau2**2) + vel[:,0,:] = (vel0*angCos).sum(axis = 1) + vel[:,1,:] = (vel0*angSin).sum(axis = 1) - pairslist1 = list(pairslist) - pairslist1.append((0,1)) - pairslist1.append((3,4)) + ind = numpy.where(numpy.isinf(vel)) + vel[ind] = numpy.nan + + return vel + + def __getPairsAutoCorr(self, pairsList, nChannels): - listMeteors1 = [] - listPowerSeries = [] - listVoltageSeries = [] - #volts has the war data + 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) - if frequency == 30e6: - timeLag = 45*10**-3 - else: - timeLag = 15*10**-3 - lag = numpy.ceil(timeLag/timeInterval) + pairsCrossCorr = range(len(pairsList)) + pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) - for i in range(len(listMeteors)): - - ###################### 3.6 - 3.7 PARAMETERS REESTIMATION ######################### - meteorAux = numpy.zeros(16) - - #Loading meteor Data (mHeight, mStart, mPeak, mEnd) - mHeight = listMeteors[i][0] - mStart = listMeteors[i][1] - mPeak = listMeteors[i][2] - mEnd = listMeteors[i][3] - - #get the volt data between the start and end times of the meteor - meteorVolts = volts[:,mStart:mEnd+1,mHeight] - meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) - - #3.6. Phase Difference estimation - phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist) - - #3.7. Phase difference removal & meteor start, peak and end times reestimated - #meteorVolts0.- all Channels, all Profiles - meteorVolts0 = volts[:,:,mHeight] - meteorThresh = noise[:,mHeight]*thresholdNoise - meteorNoise = noise[:,mHeight] - meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting - powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power - - #Times reestimation - mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0] - if mStart1.size > 0: - mStart1 = mStart1[-1] + 1 - - else: - mStart1 = mPeak - - mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1 - mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0] - if mEndDecayTime1.size == 0: - mEndDecayTime1 = powerNet0.size - else: - mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1 -# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax() - - #meteorVolts1.- all Channels, from start to end - meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1] - meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1] - if meteorVolts2.shape[1] == 0: - meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1] - meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1) - meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1) - ##################### END PARAMETERS REESTIMATION ######################### - - ##################### 3.8 PHASE DIFFERENCE REESTIMATION ######################## -# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis - if meteorVolts2.shape[1] > 0: - #Phase Difference re-estimation - phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation -# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist) - meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1]) - phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1)) - meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting - - #Phase Difference RMS - phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1))) - powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0) - #Data from Meteor - mPeak1 = powerNet1.argmax() + mStart1 - mPeakPower1 = powerNet1.max() - noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight]) - mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux - Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]) - Meteor1 = numpy.hstack((Meteor1,phaseDiffint)) - PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1] - #Vectorize - meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1] - meteorAux[7:11] = phaseDiffint[0:4] - - #Rejection Criterions - if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation - meteorAux[-1] = 17 - elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB - meteorAux[-1] = 1 - - - else: - meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd] - meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis - PowerSeries = 0 - - listMeteors1.append(meteorAux) - listPowerSeries.append(PowerSeries) - listVoltageSeries.append(meteorVolts1) - - return listMeteors1, listPowerSeries, listVoltageSeries - - def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency): + return pairsAutoCorr, pairsCrossCorr + + def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): + """ + Function that implements Spaced Antenna (SA) technique. - threshError = 10 - #Depending if it is 30 or 50 MHz - if frequency == 30e6: - timeLag = 45*10**-3 - else: - timeLag = 15*10**-3 - lag = numpy.ceil(timeLag/timeInterval) + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR - listMeteors1 = [] + Output: Winds estimation (Zonal, Meridional and Vertical) - for i in range(len(listMeteors)): - meteorPower = listPower[i] - meteorAux = listMeteors[i] - - if meteorAux[-1] == 0: - - try: - indmax = meteorPower.argmax() - indlag = indmax + lag - - y = meteorPower[indlag:] - x = numpy.arange(0, y.size)*timeLag - - #first guess - a = y[0] - tau = timeLag - #exponential fit - popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau]) - y1 = self.__exponential_function(x, *popt) - #error estimation - error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size)) - - decayTime = popt[1] - riseTime = indmax*timeInterval - meteorAux[11:13] = [decayTime, error] - - #Table items 7, 8 and 11 - if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s - meteorAux[-1] = 7 - elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time - meteorAux[-1] = 8 - if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time - meteorAux[-1] = 11 - - - except: - meteorAux[-1] = 11 - - - listMeteors1.append(meteorAux) + Parameters affected: Winds + """ + #Cross Correlation pairs obtained + pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) + pairsArray = numpy.array(pairsList)[pairsCrossCorr] + pairsSelArray = numpy.array(pairsSelected) + pairs = [] - return listMeteors1 - - #Exponential Function + #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,:] + + #--------------------------------------------------------------------- + #Metodo Directo + distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth) + winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) + 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) - def __exponential_function(self, x, a, tau): - y = a*numpy.exp(-x/tau) - return y + #--------------------------------------------------------------------- + winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) + winds = correctFactor*winds + return winds - def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval): + def __checkTime(self, currentTime, paramInterval, outputInterval): - pairslist1 = list(pairslist) - pairslist1.append((0,1)) - pairslist1.append((3,4)) - numPairs = len(pairslist1) - #Time Lag - timeLag = 45*10**-3 - c = 3e8 - lag = numpy.ceil(timeLag/timeInterval) - freq = 30e6 + dataTime = currentTime + paramInterval + deltaTime = dataTime - self.__initime - listMeteors1 = [] + if deltaTime >= outputInterval or deltaTime < 0: + self.__dataReady = True + return + + def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): + ''' + Function that implements winds estimation technique with detected meteors. - for i in range(len(listMeteors)): - meteorAux = listMeteors[i] - if meteorAux[-1] == 0: - mStart = listMeteors[i][1] - mPeak = listMeteors[i][2] - mLag = mPeak - mStart + lag - - #get the volt data between the start and end times of the meteor - meteorVolts = listVolts[i] - meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) - - #Get CCF - allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2]) - - #Method 2 - slopes = numpy.zeros(numPairs) - time = numpy.array([-2,-1,1,2])*timeInterval - angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0]) - - #Correct phases - derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1] - indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) + Input: Detected meteors, Minimum meteor quantity to wind estimation + + Output: Winds estimation (Zonal and Meridional) + + Parameters affected: Winds + ''' +# print arrayMeteor.shape + #Settings + nInt = (heightMax - heightMin)/2 +# print nInt + nInt = int(nInt) +# print nInt + winds = numpy.zeros((2,nInt))*numpy.nan + + #Filter errors + error = numpy.where(arrayMeteor[:,-1] == 0)[0] + finalMeteor = arrayMeteor[error,:] + + #Meteor Histogram + finalHeights = finalMeteor[:,2] + hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) + nMeteorsPerI = hist[0] + heightPerI = hist[1] + + #Sort of meteors + indSort = finalHeights.argsort() + finalMeteor2 = finalMeteor[indSort,:] + + # Calculating winds + ind1 = 0 + ind2 = 0 + + for i in range(nInt): + nMet = nMeteorsPerI[i] + ind1 = ind2 + ind2 = ind1 + nMet + + meteorAux = finalMeteor2[ind1:ind2,:] + + if meteorAux.shape[0] >= meteorThresh: + vel = meteorAux[:, 6] + zen = meteorAux[:, 4]*numpy.pi/180 + azim = meteorAux[:, 3]*numpy.pi/180 - if indDer[0].shape[0] > 0: - for i in range(indDer[0].shape[0]): - signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]]) - angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi - -# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]])) - for j in range(numPairs): - fit = stats.linregress(time, angAllCCF[j,:]) - slopes[j] = fit[0] + n = numpy.cos(zen) + # m = (1 - n**2)/(1 - numpy.tan(azim)**2) + # l = m*numpy.tan(azim) + l = numpy.sin(zen)*numpy.sin(azim) + m = numpy.sin(zen)*numpy.cos(azim) - #Remove Outlier -# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) -# slopes = numpy.delete(slopes,indOut) -# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) -# slopes = numpy.delete(slopes,indOut) - - radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq) - radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq) - meteorAux[-2] = radialError - meteorAux[-3] = radialVelocity + A = numpy.vstack((l, m)).transpose() + A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) + windsAux = numpy.dot(A1, vel) - #Setting Error - #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s - if numpy.abs(radialVelocity) > 200: - meteorAux[-1] = 15 - #Number 12: Poor fit to CCF variation for estimation of radial drift velocity - elif radialError > radialStdThresh: - meteorAux[-1] = 12 + winds[0,i] = windsAux[0] + winds[1,i] = windsAux[1] - listMeteors1.append(meteorAux) - return listMeteors1 + return winds, heightPerI[:-1] - def __setNewArrays(self, listMeteors, date, heiRang): - - #New arrays - arrayMeteors = numpy.array(listMeteors) - arrayParameters = numpy.zeros((len(listMeteors), 13)) + def techniqueNSM_SA(self, **kwargs): + metArray = kwargs['metArray'] + heightList = kwargs['heightList'] + timeList = kwargs['timeList'] - #Date inclusion -# date = re.findall(r'\((.*?)\)', date) -# date = date[0].split(',') -# date = map(int, date) -# -# if len(date)<6: -# date.append(0) -# -# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] -# arrayDate = numpy.tile(date, (len(listMeteors), 1)) - arrayDate = numpy.tile(date, (len(listMeteors))) + rx_location = kwargs['rx_location'] + groupList = kwargs['groupList'] + azimuth = kwargs['azimuth'] + dfactor = kwargs['dfactor'] + k = kwargs['k'] - #Meteor array -# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] -# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) + azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) + d = dist*dfactor + #Phase calculation + metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) - #Parameters Array - arrayParameters[:,0] = arrayDate #Date - arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range - arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error - arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases - arrayParameters[:,-1] = arrayMeteors[:,-1] #Error - + metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities - return arrayParameters - -# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): -# -# arrayAOA = numpy.zeros((phases.shape[0],3)) -# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList) -# -# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth) -# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1) -# arrayAOA[:,2] = cosDirError -# -# azimuthAngle = arrayAOA[:,0] -# zenithAngle = arrayAOA[:,1] -# -# #Setting Error -# #Number 3: AOA not fesible -# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0] -# error[indInvalid] = 3 -# #Number 4: Large difference in AOAs obtained from different antenna baselines -# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0] -# error[indInvalid] = 4 -# return arrayAOA, error -# -# def __getDirectionCosines(self, arrayPhase, pairsList): -# -# #Initializing some variables -# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi -# ang_aux = ang_aux.reshape(1,ang_aux.size) -# -# cosdir = numpy.zeros((arrayPhase.shape[0],2)) -# cosdir0 = numpy.zeros((arrayPhase.shape[0],2)) -# -# -# for i in range(2): -# #First Estimation -# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]] -# #Dealias -# indcsi = numpy.where(phi0_aux > numpy.pi) -# phi0_aux[indcsi] -= 2*numpy.pi -# indcsi = numpy.where(phi0_aux < -numpy.pi) -# phi0_aux[indcsi] += 2*numpy.pi -# #Direction Cosine 0 -# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5) -# -# #Most-Accurate Second Estimation -# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]] -# phi1_aux = phi1_aux.reshape(phi1_aux.size,1) -# #Direction Cosine 1 -# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5) -# -# #Searching the correct Direction Cosine -# cosdir0_aux = cosdir0[:,i] -# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) -# #Minimum Distance -# cosDiff = (cosdir1 - cosdir0_aux)**2 -# indcos = cosDiff.argmin(axis = 1) -# #Saving Value obtained -# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos] -# -# return cosdir0, cosdir -# -# def __calculateAOA(self, cosdir, azimuth): -# cosdirX = cosdir[:,0] -# cosdirY = cosdir[:,1] -# -# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi -# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east -# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose() -# -# return angles -# -# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight): -# -# Ramb = 375 #Ramb = c/(2*PRF) -# Re = 6371 #Earth Radius -# heights = numpy.zeros(Ranges.shape) -# -# R_aux = numpy.array([0,1,2])*Ramb -# R_aux = R_aux.reshape(1,R_aux.size) -# -# Ranges = Ranges.reshape(Ranges.size,1) -# -# Ri = Ranges + R_aux -# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re -# -# #Check if there is a height between 70 and 110 km -# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1) -# ind_h = numpy.where(h_bool == 1)[0] -# -# hCorr = hi[ind_h, :] -# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight)) -# -# hCorr = hi[ind_hCorr] -# heights[ind_h] = hCorr -# -# #Setting Error -# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km -# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km -# -# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] -# error[indInvalid2] = 14 -# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0] -# error[indInvalid1] = 13 -# -# return heights, error + velEst = numpy.zeros((heightList.size,2))*numpy.nan + azimuth1 = azimuth1*numpy.pi/180 + + for i in range(heightList.size): + h = heightList[i] + indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] + metHeight = metArray1[indH,:] + if metHeight.shape[0] >= 2: + velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities + iazim = metHeight[:,1].astype(int) + azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths + A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux))) + A = numpy.asmatrix(A) + A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() + velHor = numpy.dot(A1,velAux) + + velEst[i,:] = numpy.squeeze(velHor) + return velEst - def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None): + def __getPhaseSlope(self, metArray, heightList, timeList): + meteorList = [] + #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 + #Putting back together the meteor matrix + utctime = metArray[:,0] + uniqueTime = numpy.unique(utctime) - ''' - Function GetMoments() + phaseDerThresh = 0.5 + ippSeconds = timeList[1] - timeList[0] + sec = numpy.where(timeList>1)[0][0] + nPairs = metArray.shape[1] - 6 + nHeights = len(heightList) - Input: - Output: - Variables modified: - ''' - if path != None: - sys.path.append(path) - self.dataOut.library = importlib.import_module(file) + for t in uniqueTime: + metArray1 = metArray[utctime==t,:] +# phaseDerThresh = numpy.pi/4 #reducir Phase thresh + tmet = metArray1[:,1].astype(int) + hmet = metArray1[:,2].astype(int) + + metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) + metPhase[:,:] = numpy.nan + metPhase[:,hmet,tmet] = metArray1[:,6:].T + + #Delete short trails + metBool = ~numpy.isnan(metPhase[0,:,:]) + heightVect = numpy.sum(metBool, axis = 1) + metBool[heightVect phaseDerThresh)) + metPhase[phDerAux] = numpy.nan + + #--------------------------METEOR DETECTION ----------------------------------------- + indMet = numpy.where(numpy.any(metBool,axis=1))[0] + + for p in numpy.arange(nPairs): + phase = metPhase[p,:,:] + phDer = metDer[p,:,:] + + for h in indMet: + height = heightList[h] + phase1 = phase[h,:] #82 + phDer1 = phDer[h,:] + + phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap + + indValid = numpy.where(~numpy.isnan(phase1))[0] + initMet = indValid[0] + endMet = 0 + + for i in range(len(indValid)-1): + + #Time difference + inow = indValid[i] + inext = indValid[i+1] + idiff = inext - inow + #Phase difference + phDiff = numpy.abs(phase1[inext] - phase1[inow]) + + if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor + sizeTrail = inow - initMet + 1 + if sizeTrail>3*sec: #Too short meteors + x = numpy.arange(initMet,inow+1)*ippSeconds + y = phase1[initMet:inow+1] + ynnan = ~numpy.isnan(y) + x = x[ynnan] + y = y[ynnan] + slope, intercept, r_value, p_value, std_err = stats.linregress(x,y) + ylin = x*slope + intercept + rsq = r_value**2 + if rsq > 0.5: + vel = slope#*height*1000/(k*d) + estAux = numpy.array([utctime,p,height, vel, rsq]) + meteorList.append(estAux) + initMet = inext + metArray2 = numpy.array(meteorList) - #To be inserted as a parameter - groupArray = numpy.array(groupList) -# groupArray = numpy.array([[0,1],[2,3]]) - self.dataOut.groupList = groupArray + return metArray2 + + def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): - nGroups = groupArray.shape[0] - nChannels = self.dataIn.nChannels - nHeights=self.dataIn.heightList.size + azimuth1 = numpy.zeros(len(pairslist)) + dist = numpy.zeros(len(pairslist)) - #Parameters Array - self.dataOut.data_param = None + for i in range(len(rx_location)): + ch0 = pairslist[i][0] + ch1 = pairslist[i][1] + + diffX = rx_location[ch0][0] - rx_location[ch1][0] + diffY = rx_location[ch0][1] - rx_location[ch1][1] + azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi + dist[i] = numpy.sqrt(diffX**2 + diffY**2) - #Set constants - constants = self.dataOut.library.setConstants(self.dataIn) - self.dataOut.constants = constants - M = self.dataIn.normFactor - N = self.dataIn.nFFTPoints - ippSeconds = self.dataIn.ippSeconds - K = self.dataIn.nIncohInt - pairsArray = numpy.array(self.dataIn.pairsList) + azimuth1 -= azimuth0 + return azimuth1, dist + + def techniqueNSM_DBS(self, **kwargs): + metArray = kwargs['metArray'] + heightList = kwargs['heightList'] + timeList = kwargs['timeList'] + zenithList = kwargs['zenithList'] + nChan = numpy.max(cmet) + 1 + nHeights = len(heightList) - #List of possible combinations - listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) - indCross = numpy.zeros(len(list(listComb)), dtype = 'int') + utctime = metArray[:,0] + cmet = metArray[:,1] + hmet = metArray1[:,3].astype(int) + h1met = heightList[hmet]*zenithList[cmet] + vmet = metArray1[:,5] - if getSNR: - listChannels = groupArray.reshape((groupArray.size)) - listChannels.sort() - noise = self.dataIn.getNoise() - self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) + for i in range(nHeights - 1): + hmin = heightList[i] + hmax = heightList[i + 1] + + vthisH = vmet[(h1met>=hmin) & (h1met0)and(error1[3]<5): - p0 = self.dataOut.data_param[i,:,h-1] - else: - p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i)) + dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) + dataOut.flagNoData = False + self.__buffer = None - try: - #Least Squares - minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) -# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) - #Chi square error - error0 = numpy.sum(infodict['fvec']**2)/(2*N) - #Error with Jacobian - error1 = self.dataOut.library.errorFunction(minp,constants,LT) - except: - minp = p0*numpy.nan - error0 = numpy.nan - error1 = p0*numpy.nan - - #Save - if self.dataOut.data_param == None: - self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan - self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan + elif technique == 'Meteors1': + dataOut.flagNoData = True + self.__dataReady = False + + if kwargs.has_key('nMins'): + nMins = kwargs['nMins'] + else: nMins = 20 + if kwargs.has_key('rx_location'): + rx_location = kwargs['rx_location'] + else: rx_location = [(0,1),(1,1),(1,0)] + if kwargs.has_key('azimuth'): + azimuth = kwargs['azimuth'] + else: azimuth = 51 + if kwargs.has_key('dfactor'): + dfactor = kwargs['dfactor'] + if kwargs.has_key('mode'): + mode = kwargs['mode'] + else: mode = 'SA' + + #Borrar luego esto + if dataOut.groupList == None: + dataOut.groupList = [(0,1),(0,2),(1,2)] + groupList = dataOut.groupList + C = 3e8 + freq = 50e6 + lamb = C/freq + k = 2*numpy.pi/lamb + + timeList = dataOut.abscissaList + heightList = dataOut.heightList + + if self.__isConfig == False: + dataOut.outputInterval = nMins*60 +# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) + #Get Initial LTC time + initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) + minuteAux = initime.minute + minuteNew = int(numpy.floor(minuteAux/nMins)*nMins) + self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + + self.__isConfig = True - self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) - self.dataOut.data_param[i,:,h] = minp + if self.__buffer == None: + self.__buffer = dataOut.data_param + self.__firstdata = copy.copy(dataOut) + + else: + self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + + self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready + + if self.__dataReady: + dataOut.utctimeInit = self.__initime + self.__initime += dataOut.outputInterval #to erase time offset + + metArray = self.__buffer + if mode == 'SA': + dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) + elif mode == 'DBS': + dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList) + dataOut.data_output = dataOut.data_output.T + dataOut.flagNoData = False + self.__buffer = None + + return + +class EWDriftsEstimation(Operation): + + def __init__(self): + Operation.__init__(self) + + def __correctValues(self, heiRang, phi, velRadial, SNR): + listPhi = phi.tolist() + maxid = listPhi.index(max(listPhi)) + minid = listPhi.index(min(listPhi)) + + rango = range(len(phi)) + # rango = numpy.delete(rango,maxid) + + heiRang1 = heiRang*math.cos(phi[maxid]) + heiRangAux = heiRang*math.cos(phi[minid]) + indOut = (heiRang1 < heiRangAux[0]).nonzero() + heiRang1 = numpy.delete(heiRang1,indOut) + + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) + SNR1 = numpy.zeros([len(phi),len(heiRang1)]) + + for i in rango: + x = heiRang*math.cos(phi[i]) + y1 = velRadial[i,:] + f1 = interpolate.interp1d(x,y1,kind = 'cubic') + + x1 = heiRang1 + y11 = f1(x1) + + y2 = SNR[i,:] + f2 = interpolate.interp1d(x,y2,kind = 'cubic') + y21 = f2(x1) + + velRadial1[i,:] = y11 + SNR1[i,:] = y21 + + return heiRang1, velRadial1, SNR1 + + def run(self, dataOut, zenith, zenithCorrection): + heiRang = dataOut.heightList + velRadial = dataOut.data_param[:,3,:] + SNR = dataOut.data_SNR + + zenith = numpy.array(zenith) + zenith -= zenithCorrection + zenith *= numpy.pi/180 + + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) + + alp = zenith[0] + bet = zenith[1] + + w_w = velRadial1[0,:] + w_e = velRadial1[1,:] + + w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) + u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) + + winds = numpy.vstack((u,w)) + + dataOut.heightList = heiRang1 + dataOut.data_output = winds + dataOut.data_SNR = SNR1 + + dataOut.utctimeInit = dataOut.utctime + dataOut.outputInterval = dataOut.timeInterval return - - def __residFunction(self, p, dp, LT, constants): - fm = self.dataOut.library.modelFunction(p, constants) - fmp=numpy.dot(LT,fm) - - return dp-fmp +#--------------- Non Specular Meteor ---------------- - def __getSNR(self, z, noise): - - avg = numpy.average(z, axis=1) - SNR = (avg.T-noise)/noise - SNR = SNR.T - return SNR - - def __chisq(p,chindex,hindex): - #similar to Resid but calculates CHI**2 - [LT,d,fm]=setupLTdfm(p,chindex,hindex) - dp=numpy.dot(LT,d) - fmp=numpy.dot(LT,fm) - chisq=numpy.dot((dp-fmp).T,(dp-fmp)) - return chisq - - def NonSpecularMeteorDetection(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False): +class NonSpecularMeteorDetection(Operation): + + def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False): data_acf = self.dataOut.data_pre[0] data_ccf = self.dataOut.data_pre[1] @@ -1485,806 +1395,802 @@ class ParametersProc(ProcessingUnit): coordMet = numpy.where(boolMetFin) cmet = coordMet[0] - tmet = coordMet[1] - hmet = coordMet[2] - - data_param = numpy.zeros((tmet.size, 7)) - data_param[:,0] = utctime - data_param[:,1] = cmet - data_param[:,2] = tmet - data_param[:,3] = hmet - data_param[:,4] = SNR[cmet,tmet,hmet].T - data_param[:,5] = velRad[cmet,tmet,hmet].T - data_param[:,6] = spcWidth[cmet,tmet,hmet].T - -# self.dataOut.data_param = data_int - if len(data_param) == 0: - self.dataOut.flagNoData = True - else: - self.dataOut.data_param = data_param - - def __erase_small(self, binArray, threshX, threshY): - labarray, numfeat = ndimage.measurements.label(binArray) - binArray1 = numpy.copy(binArray) - - for i in range(1,numfeat + 1): - auxBin = (labarray==i) - auxSize = auxBin.sum() - - x,y = numpy.where(auxBin) - widthX = x.max() - x.min() - widthY = y.max() - y.min() - - #width X: 3 seg -> 12.5*3 - #width Y: - - if (auxSize < 50) or (widthX < threshX) or (widthY < threshY): - binArray1[auxBin] = False - - return binArray1 - - def WeirdEcho(self): -# data_pre = self.dataOut.data_pre -# nHeights = self.dataOut.nHeights -# # nProfiles = self.dataOut.data_pre.shape[1] -# # data_param = numpy.zeros((len(pairslist),nProfiles,nHeights)) -# data_param = numpy.zeros((len(pairslist),nHeights)) -# -# for i in range(len(pairslist)): -# chan0 = data_pre[pairslist[i][0],:] -# chan1 = data_pre[pairslist[i][1],:] -# #calcular correlacion cruzada -# #magnitud es coherencia -# #fase es dif fase -# correl = chan0*numpy.conj(chan1) -# coherence = numpy.abs(correl)/(numpy.abs(chan0)*numpy.abs(chan1)) -# phase = numpy.angle(correl) -# # data_param[2*i,:,:] = coherence -# data_param[i,:] = phase -# -# self.dataOut.data_param = data_param -# self.dataOut.groupList = pairslist - data_cspc = self.dataOut.data_pre[1] - ccf = numpy.average(data_cspc,axis=1) - phases = numpy.angle(ccf).T - - meteorOps = MeteorOperations() - pairsList = ((0,1),(2,3)) - jph = numpy.array([0,0,0,0]) - AOAthresh = numpy.pi/8 - azimuth = 45 - error = numpy.zeros((phases.shape[0],1)) - AOA,error = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) - self.dataOut.data_param = AOA.T - -class WindProfiler(Operation): - - __isConfig = False - - __initime = None - __lastdatatime = None - __integrationtime = None - - __buffer = None - - __dataReady = False - - __firstdata = None - - n = None - - def __init__(self): - Operation.__init__(self) - - def __calculateCosDir(self, elev, azim): - zen = (90 - elev)*numpy.pi/180 - azim = azim*numpy.pi/180 - cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2))) - cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2) - - signX = numpy.sign(numpy.cos(azim)) - signY = numpy.sign(numpy.sin(azim)) - - cosDirX = numpy.copysign(cosDirX, signX) - cosDirY = numpy.copysign(cosDirY, signY) - return cosDirX, cosDirY - - def __calculateAngles(self, theta_x, theta_y, azimuth): - - dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) - zenith_arr = numpy.arccos(dir_cosw) - azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 - - dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) - dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) - - return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw - - def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): - -# - if horOnly: - A = numpy.c_[dir_cosu,dir_cosv] - else: - A = numpy.c_[dir_cosu,dir_cosv,dir_cosw] - A = numpy.asmatrix(A) - A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose() - - return A1 - - def __correctValues(self, heiRang, phi, velRadial, SNR): - listPhi = phi.tolist() - maxid = listPhi.index(max(listPhi)) - minid = listPhi.index(min(listPhi)) - - rango = range(len(phi)) - # rango = numpy.delete(rango,maxid) - - heiRang1 = heiRang*math.cos(phi[maxid]) - heiRangAux = heiRang*math.cos(phi[minid]) - indOut = (heiRang1 < heiRangAux[0]).nonzero() - heiRang1 = numpy.delete(heiRang1,indOut) - - velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) - SNR1 = numpy.zeros([len(phi),len(heiRang1)]) - - for i in rango: - x = heiRang*math.cos(phi[i]) - y1 = velRadial[i,:] - f1 = interpolate.interp1d(x,y1,kind = 'cubic') - - x1 = heiRang1 - y11 = f1(x1) - - y2 = SNR[i,:] - f2 = interpolate.interp1d(x,y2,kind = 'cubic') - y21 = f2(x1) - - velRadial1[i,:] = y11 - SNR1[i,:] = y21 - - return heiRang1, velRadial1, SNR1 - - 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): - """ - Function that implements Doppler Beam Swinging (DBS) technique. - - Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, - Direction correction (if necessary), Ranges and SNR - - Output: Winds estimation (Zonal, Meridional and Vertical) - - Parameters affected: Winds, height range, SNR - """ - azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth) - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0) - A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) - - #Calculo de Componentes de la velocidad con DBS - winds = self.__calculateVelUVW(A,velRadial1) - - return winds, heiRang1, SNR1 - - def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None): - - posx = numpy.asarray(posx) - posy = numpy.asarray(posy) - - #Rotacion Inversa para alinear con el azimuth - if azimuth!= None: - azimuth = azimuth*math.pi/180 - posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth) - posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth) - else: - posx1 = posx - posy1 = posy - - #Calculo de Distancias - distx = numpy.zeros(pairsCrossCorr.size) - disty = numpy.zeros(pairsCrossCorr.size) - dist = numpy.zeros(pairsCrossCorr.size) - ang = numpy.zeros(pairsCrossCorr.size) - - for i in range(pairsCrossCorr.size): - distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]] - disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]] - dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) - ang[i] = numpy.arctan2(disty[i],distx[i]) - #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]] + tmet = coordMet[1] + hmet = coordMet[2] - return distx,disty, dist1,ang1 - - def __calculateVelVer(self, phase, lagTRange, _lambda): + data_param = numpy.zeros((tmet.size, 7)) + data_param[:,0] = utctime + data_param[:,1] = cmet + data_param[:,2] = tmet + data_param[:,3] = hmet + data_param[:,4] = SNR[cmet,tmet,hmet].T + data_param[:,5] = velRad[cmet,tmet,hmet].T + data_param[:,6] = spcWidth[cmet,tmet,hmet].T + +# self.dataOut.data_param = data_int + if len(data_param) == 0: + self.dataOut.flagNoData = True + else: + self.dataOut.data_param = data_param - Ts = lagTRange[1] - lagTRange[0] - velW = -_lambda*phase/(4*math.pi*Ts) - - return velW - - def __calculateVelHorDir(self, dist, tau1, tau2, ang): - nPairs = tau1.shape[0] - vel = numpy.zeros((nPairs,3,tau1.shape[2])) + def __erase_small(self, binArray, threshX, threshY): + labarray, numfeat = ndimage.measurements.label(binArray) + binArray1 = numpy.copy(binArray) - angCos = numpy.cos(ang) - angSin = numpy.sin(ang) + for i in range(1,numfeat + 1): + auxBin = (labarray==i) + auxSize = auxBin.sum() + + x,y = numpy.where(auxBin) + widthX = x.max() - x.min() + widthY = y.max() - y.min() + + #width X: 3 seg -> 12.5*3 + #width Y: + + if (auxSize < 50) or (widthX < threshX) or (widthY < threshY): + binArray1[auxBin] = False + + return binArray1 + +#--------------- Specular Meteor ---------------- + +class MeteorDetection(Operation): + ''' + Function DetectMeteors() + Project developed with paper: + HOLDSWORTH ET AL. 2004 + + Input: + self.dataOut.data_pre - vel0 = dist*tau1/(2*tau2**2) - vel[:,0,:] = (vel0*angCos).sum(axis = 1) - vel[:,1,:] = (vel0*angSin).sum(axis = 1) + centerReceiverIndex: From the channels, which is the center receiver + + hei_ref: Height reference for the Beacon signal extraction + tauindex: + predefinedPhaseShifts: Predefined phase offset for the voltge signals + + cohDetection: Whether to user Coherent detection or not + cohDet_timeStep: Coherent Detection calculation time step + cohDet_thresh: Coherent Detection phase threshold to correct phases + + noise_timeStep: Noise calculation time step + noise_multiple: Noise multiple to define signal threshold + + multDet_timeLimit: Multiple Detection Removal time limit in seconds + multDet_rangeLimit: Multiple Detection Removal range limit in km + + phaseThresh: Maximum phase difference between receiver to be consider a meteor + SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor + + hmin: Minimum Height of the meteor to use it in the further wind estimations + hmax: Maximum Height of the meteor to use it in the further wind estimations + azimuth: Azimuth angle correction + + Affected: + self.dataOut.data_param - ind = numpy.where(numpy.isinf(vel)) - vel[ind] = numpy.nan - - return vel - - def __getPairsAutoCorr(self, pairsList, nChannels): - - pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan + Rejection Criteria (Errors): + 0: No error; analysis OK + 1: SNR < SNR threshold + 2: angle of arrival (AOA) ambiguously determined + 3: AOA estimate not feasible + 4: Large difference in AOAs obtained from different antenna baselines + 5: echo at start or end of time series + 6: echo less than 5 examples long; too short for analysis + 7: echo rise exceeds 0.3s + 8: echo decay time less than twice rise time + 9: large power level before echo + 10: large power level after echo + 11: poor fit to amplitude for estimation of decay time + 12: poor fit to CCF phase variation for estimation of radial drift velocity + 13: height unresolvable echo: not valid height within 70 to 110 km + 14: height ambiguous echo: more then one possible height within 70 to 110 km + 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s + 16: oscilatory echo, indicating event most likely not an underdense echo - 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) + 17: phase difference in meteor Reestimation - pairsCrossCorr = range(len(pairsList)) - pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) + Data Storage: + Meteors for Wind Estimation (8): + Utc Time | Range Height + Azimuth Zenith errorCosDir + VelRad errorVelRad + Phase0 Phase1 Phase2 Phase3 + TypeError - return pairsAutoCorr, pairsCrossCorr + ''' - def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): - """ - Function that implements Spaced Antenna (SA) technique. - - Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, - Direction correction (if necessary), Ranges and SNR + def run(self, dataOut, hei_ref = None, tauindex = 0, + phaseOffsets = None, + cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, + noise_timeStep = 4, noise_multiple = 4, + multDet_timeLimit = 1, multDet_rangeLimit = 3, + phaseThresh = 20, SNRThresh = 5, + hmin = 50, hmax=150, azimuth = 0, + channelPositions = None) : - Output: Winds estimation (Zonal, Meridional and Vertical) + + #Getting Pairslist + if channelPositions == None: +# channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T + channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella + meteorOps = MeteorOperations() + pairslist0, distances = meteorOps.getPhasePairs(channelPositions) - Parameters affected: Winds - """ - #Cross Correlation pairs obtained - pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) - pairsArray = numpy.array(pairsList)[pairsCrossCorr] - pairsSelArray = numpy.array(pairsSelected) - pairs = [] + #Get Beacon signal + newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex]) - #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)) + if hei_ref != None: + newheis = numpy.where(self.dataOut.heightList>hei_ref) - indtau = tau.shape[0]/2 - tau1 = tau[:indtau,:] - tau2 = tau[indtau:-1,:] - tau1 = tau1[pairs,:] - tau2 = tau2[pairs,:] - phase1 = tau[-1,:] + heiRang = dataOut.getHeiRange() + + #****************REMOVING HARDWARE PHASE DIFFERENCES*************** + # see if the user put in pre defined phase shifts + voltsPShift = self.dataOut.data_pre.copy() - #--------------------------------------------------------------------- - #Metodo Directo - distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth) - winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) - 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) +# if predefinedPhaseShifts != None: +# hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 +# +# # elif beaconPhaseShifts: +# # #get hardware phase shifts using beacon signal +# # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) +# # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) +# +# else: +# hardwarePhaseShifts = numpy.zeros(5) +# +# voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex') +# for i in range(self.dataOut.data_pre.shape[0]): +# voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i]) - #--------------------------------------------------------------------- - winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) - winds = correctFactor*winds - return winds - - def __checkTime(self, currentTime, paramInterval, outputInterval): - - dataTime = currentTime + paramInterval - deltaTime = dataTime - self.__initime - - if deltaTime >= outputInterval or deltaTime < 0: - self.__dataReady = True - return + #******************END OF REMOVING HARDWARE PHASE DIFFERENCES********* - def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): - ''' - Function that implements winds estimation technique with detected meteors. - - Input: Detected meteors, Minimum meteor quantity to wind estimation - - Output: Winds estimation (Zonal and Meridional) - - Parameters affected: Winds - ''' -# print arrayMeteor.shape - #Settings - nInt = (heightMax - heightMin)/2 -# print nInt - nInt = int(nInt) -# print nInt - winds = numpy.zeros((2,nInt))*numpy.nan - - #Filter errors - error = numpy.where(arrayMeteor[:,-1] == 0)[0] - finalMeteor = arrayMeteor[error,:] + #Remove DC + voltsDC = numpy.mean(voltsPShift,1) + voltsDC = numpy.mean(voltsDC,1) + for i in range(voltsDC.shape[0]): + voltsPShift[i] = voltsPShift[i] - voltsDC[i] + + #Don't considerate last heights, theyre used to calculate Hardware Phase Shift + voltsPShift = voltsPShift[:,:,:newheis[0][0]] - #Meteor Histogram - finalHeights = finalMeteor[:,2] - hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) - nMeteorsPerI = hist[0] - heightPerI = hist[1] + #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) ********** + #Coherent Detection + if cohDetection: + #use coherent detection to get the net power + cohDet_thresh = cohDet_thresh*numpy.pi/180 + voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist0, cohDet_thresh) - #Sort of meteors - indSort = finalHeights.argsort() - finalMeteor2 = finalMeteor[indSort,:] + #Non-coherent detection! + powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0) + #********** END OF COH/NON-COH POWER CALCULATION********************** + + #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS **************** + #Get noise + noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval) +# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval) + #Get signal threshold + signalThresh = noise_multiple*noise + #Meteor echoes detection + listMeteors = self.__findMeteors(powerNet, signalThresh) + #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********** - # Calculating winds - ind1 = 0 - ind2 = 0 + #************** REMOVE MULTIPLE DETECTIONS (3.5) *************************** + #Parameters + heiRange = self.dataOut.getHeiRange() + rangeInterval = heiRange[1] - heiRange[0] + rangeLimit = multDet_rangeLimit/rangeInterval + timeLimit = multDet_timeLimit/self.dataOut.timeInterval + #Multiple detection removals + listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit) + #************ END OF REMOVE MULTIPLE DETECTIONS ********************** - for i in range(nInt): - nMet = nMeteorsPerI[i] - ind1 = ind2 - ind2 = ind1 + nMet + #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ******************** + #Parameters + phaseThresh = phaseThresh*numpy.pi/180 + thresh = [phaseThresh, noise_multiple, SNRThresh] + #Meteor reestimation (Errors N 1, 6, 12, 17) + listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency) +# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise) + #Estimation of decay times (Errors N 7, 8, 11) + listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency) + #******************* END OF METEOR REESTIMATION ******************* + + #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) ************************** + #Calculating Radial Velocity (Error N 15) + radialStdThresh = 10 + listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, self.dataOut.timeInterval) + + if len(listMeteors4) > 0: + #Setting New Array + date = self.dataOut.utctime + arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) - meteorAux = finalMeteor2[ind1:ind2,:] + #Correcting phase offset + if phaseOffsets != None: + phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 + arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) - if meteorAux.shape[0] >= meteorThresh: - vel = meteorAux[:, 6] - zen = meteorAux[:, 4]*numpy.pi/180 - azim = meteorAux[:, 3]*numpy.pi/180 - - n = numpy.cos(zen) - # m = (1 - n**2)/(1 - numpy.tan(azim)**2) - # l = m*numpy.tan(azim) - l = numpy.sin(zen)*numpy.sin(azim) - m = numpy.sin(zen)*numpy.cos(azim) - - A = numpy.vstack((l, m)).transpose() - A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) - windsAux = numpy.dot(A1, vel) - - winds[0,i] = windsAux[0] - winds[1,i] = windsAux[1] + #Second Pairslist + pairsList = [] + pairx = (0,1) + pairy = (2,3) + pairsList.append(pairx) + pairsList.append(pairy) + + jph = numpy.array([0,0,0,0]) + h = (hmin,hmax) + arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) + +# #Calculate AOA (Error N 3, 4) +# #JONES ET AL. 1998 +# error = arrayParameters[:,-1] +# AOAthresh = numpy.pi/8 +# phases = -arrayParameters[:,9:13] +# arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) +# +# #Calculate Heights (Error N 13 and 14) +# error = arrayParameters[:,-1] +# Ranges = arrayParameters[:,2] +# zenith = arrayParameters[:,5] +# arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) +# error = arrayParameters[:,-1] + #********************* END OF PARAMETERS CALCULATION ************************** - return winds, heightPerI[:-1] - - def techniqueNSM_SA(self, **kwargs): - metArray = kwargs['metArray'] - heightList = kwargs['heightList'] - timeList = kwargs['timeList'] + #***************************+ PASS DATA TO NEXT STEP ********************** +# arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) + self.dataOut.data_param = arrayParameters + + if arrayParameters == None: + self.dataOut.flagNoData = True + else: + self.dataOut.flagNoData = True + + return - rx_location = kwargs['rx_location'] - groupList = kwargs['groupList'] - azimuth = kwargs['azimuth'] - dfactor = kwargs['dfactor'] - k = kwargs['k'] + def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n): + + minIndex = min(newheis[0]) + maxIndex = max(newheis[0]) - azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth) - d = dist*dfactor - #Phase calculation - metArray1 = self.__getPhaseSlope(metArray, heightList, timeList) + voltage = voltage0[:,:,minIndex:maxIndex+1] + nLength = voltage.shape[1]/n + nMin = 0 + nMax = 0 + phaseOffset = numpy.zeros((len(pairslist),n)) - metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities + for i in range(n): + nMax += nLength + phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0])) + phaseCCF = numpy.mean(phaseCCF, axis = 2) + phaseOffset[:,i] = phaseCCF.transpose() + nMin = nMax +# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist) - velEst = numpy.zeros((heightList.size,2))*numpy.nan - azimuth1 = azimuth1*numpy.pi/180 + #Remove Outliers + factor = 2 + wt = phaseOffset - signal.medfilt(phaseOffset,(1,5)) + dw = numpy.std(wt,axis = 1) + dw = dw.reshape((dw.size,1)) + ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) + phaseOffset[ind] = numpy.nan + phaseOffset = stats.nanmean(phaseOffset, axis=1) - for i in range(heightList.size): - h = heightList[i] - indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0] - metHeight = metArray1[indH,:] - if metHeight.shape[0] >= 2: - velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities - iazim = metHeight[:,1].astype(int) - azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths - A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux))) - A = numpy.asmatrix(A) - A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose() - velHor = numpy.dot(A1,velAux) - - velEst[i,:] = numpy.squeeze(velHor) - return velEst + return phaseOffset - def __getPhaseSlope(self, metArray, heightList, timeList): - meteorList = [] - #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2 - #Putting back together the meteor matrix - utctime = metArray[:,0] - uniqueTime = numpy.unique(utctime) - - phaseDerThresh = 0.5 - ippSeconds = timeList[1] - timeList[0] - sec = numpy.where(timeList>1)[0][0] - nPairs = metArray.shape[1] - 6 - nHeights = len(heightList) + def __shiftPhase(self, data, phaseShift): + #this will shift the phase of a complex number + dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j) + return dataShifted + + def __estimatePhaseDifference(self, array, pairslist): + nChannel = array.shape[0] + nHeights = array.shape[2] + numPairs = len(pairslist) +# phaseCCF = numpy.zeros((nChannel, 5, nHeights)) + phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2])) - for t in uniqueTime: - metArray1 = metArray[utctime==t,:] -# phaseDerThresh = numpy.pi/4 #reducir Phase thresh - tmet = metArray1[:,1].astype(int) - hmet = metArray1[:,2].astype(int) - - metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1)) - metPhase[:,:] = numpy.nan - metPhase[:,hmet,tmet] = metArray1[:,6:].T - - #Delete short trails - metBool = ~numpy.isnan(metPhase[0,:,:]) - heightVect = numpy.sum(metBool, axis = 1) - metBool[heightVect phaseDerThresh)) - metPhase[phDerAux] = numpy.nan - - #--------------------------METEOR DETECTION ----------------------------------------- - indMet = numpy.where(numpy.any(metBool,axis=1))[0] - - for p in numpy.arange(nPairs): - phase = metPhase[p,:,:] - phDer = metDer[p,:,:] - - for h in indMet: - height = heightList[h] - phase1 = phase[h,:] #82 - phDer1 = phDer[h,:] - - phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap - - indValid = numpy.where(~numpy.isnan(phase1))[0] - initMet = indValid[0] - endMet = 0 - - for i in range(len(indValid)-1): - - #Time difference - inow = indValid[i] - inext = indValid[i+1] - idiff = inext - inow - #Phase difference - phDiff = numpy.abs(phase1[inext] - phase1[inow]) - - if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor - sizeTrail = inow - initMet + 1 - if sizeTrail>3*sec: #Too short meteors - x = numpy.arange(initMet,inow+1)*ippSeconds - y = phase1[initMet:inow+1] - ynnan = ~numpy.isnan(y) - x = x[ynnan] - y = y[ynnan] - slope, intercept, r_value, p_value, std_err = stats.linregress(x,y) - ylin = x*slope + intercept - rsq = r_value**2 - if rsq > 0.5: - vel = slope#*height*1000/(k*d) - estAux = numpy.array([utctime,p,height, vel, rsq]) - meteorList.append(estAux) - initMet = inext - metArray2 = numpy.array(meteorList) + #Correct phases + derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:] + indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) - return metArray2 + if indDer[0].shape[0] > 0: + for i in range(indDer[0].shape[0]): + signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]]) + phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi - def __calculateAzimuth1(self, rx_location, pairslist, azimuth0): +# for j in range(numSides): +# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2]) +# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux) +# + #Linear + phaseInt = numpy.zeros((numPairs,1)) + angAllCCF = phaseCCF[:,[0,1,3,4],0] + for j in range(numPairs): + fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:]) + phaseInt[j] = fit[1] + #Phase Differences + phaseDiff = phaseInt - phaseCCF[:,2,:] + phaseArrival = phaseInt.reshape(phaseInt.size) - azimuth1 = numpy.zeros(len(pairslist)) - dist = numpy.zeros(len(pairslist)) + #Dealias + phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival)) +# indAlias = numpy.where(phaseArrival > numpy.pi) +# phaseArrival[indAlias] -= 2*numpy.pi +# indAlias = numpy.where(phaseArrival < -numpy.pi) +# phaseArrival[indAlias] += 2*numpy.pi - for i in range(len(rx_location)): - ch0 = pairslist[i][0] - ch1 = pairslist[i][1] + return phaseDiff, phaseArrival + + def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh): + #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power + #find the phase shifts of each channel over 1 second intervals + #only look at ranges below the beacon signal + numProfPerBlock = numpy.ceil(timeSegment/timeInterval) + numBlocks = int(volts.shape[1]/numProfPerBlock) + numHeights = volts.shape[2] + nChannel = volts.shape[0] + voltsCohDet = volts.copy() + + pairsarray = numpy.array(pairslist) + indSides = pairsarray[:,1] +# indSides = numpy.array(range(nChannel)) +# indSides = numpy.delete(indSides, indCenter) +# +# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0) + listBlocks = numpy.array_split(volts, numBlocks, 1) + + startInd = 0 + endInd = 0 - diffX = rx_location[ch0][0] - rx_location[ch1][0] - diffY = rx_location[ch0][1] - rx_location[ch1][1] - azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi - dist[i] = numpy.sqrt(diffX**2 + diffY**2) + for i in range(numBlocks): + startInd = endInd + endInd = endInd + listBlocks[i].shape[1] + + arrayBlock = listBlocks[i] +# arrayBlockCenter = listCenter[i] + + #Estimate the Phase Difference + phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist) + #Phase Difference RMS + arrayPhaseRMS = numpy.abs(phaseDiff) + phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0) + indPhase = numpy.where(phaseRMSaux==4) + #Shifting + if indPhase[0].shape[0] > 0: + for j in range(indSides.size): + arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose()) + voltsCohDet[:,startInd:endInd,:] = arrayBlock + + return voltsCohDet + + def __calculateCCF(self, volts, pairslist ,laglist): - azimuth1 -= azimuth0 - return azimuth1, dist + nHeights = volts.shape[2] + nPoints = volts.shape[1] + voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex') + + for i in range(len(pairslist)): + volts1 = volts[pairslist[i][0]] + volts2 = volts[pairslist[i][1]] + + for t in range(len(laglist)): + idxT = laglist[t] + if idxT >= 0: + vStacked = numpy.vstack((volts2[idxT:,:], + numpy.zeros((idxT, nHeights),dtype='complex'))) + else: + vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'), + volts2[:(nPoints + idxT),:])) + voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0) - def techniqueNSM_DBS(self, **kwargs): - metArray = kwargs['metArray'] - heightList = kwargs['heightList'] - timeList = kwargs['timeList'] - zenithList = kwargs['zenithList'] - nChan = numpy.max(cmet) + 1 - nHeights = len(heightList) + vStacked = None + return voltsCCF - utctime = metArray[:,0] - cmet = metArray[:,1] - hmet = metArray1[:,3].astype(int) - h1met = heightList[hmet]*zenithList[cmet] - vmet = metArray1[:,5] + def __getNoise(self, power, timeSegment, timeInterval): + numProfPerBlock = numpy.ceil(timeSegment/timeInterval) + numBlocks = int(power.shape[0]/numProfPerBlock) + numHeights = power.shape[1] + + listPower = numpy.array_split(power, numBlocks, 0) + noise = numpy.zeros((power.shape[0], power.shape[1])) + noise1 = numpy.zeros((power.shape[0], power.shape[1])) - for i in range(nHeights - 1): - hmin = heightList[i] - hmax = heightList[i + 1] + startInd = 0 + endInd = 0 + + for i in range(numBlocks): #split por canal + startInd = endInd + endInd = endInd + listPower[i].shape[0] - vthisH = vmet[(h1met>=hmin) & (h1met threshAux)[0] + indDNthresh = numpy.where(powerAux <= threshAux)[0] - return data_output + j = 0 - def run(self, dataOut, technique, **kwargs): + while (j < indUPthresh.size - 2): + if (indUPthresh[j + 2] == indUPthresh[j] + 2): + indDNAux = numpy.where(indDNthresh > indUPthresh[j]) + indDNthresh = indDNthresh[indDNAux] + + if (indDNthresh.size > 0): + indEnd = indDNthresh[0] - 1 + indInit = indUPthresh[j] + + meteor = powerAux[indInit:indEnd + 1] + indPeak = meteor.argmax() + indInit + FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0))) + + listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!! + j = numpy.where(indUPthresh == indEnd)[0] + 1 + else: j+=1 + else: j+=1 + + return listMeteors + + def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit): - param = dataOut.data_param - if dataOut.abscissaList != None: - absc = dataOut.abscissaList[:-1] - noise = dataOut.noise - heightList = dataOut.heightList - SNR = dataOut.data_SNR + arrayMeteors = numpy.asarray(listMeteors) + listMeteors1 = [] - if technique == 'DBS': + while arrayMeteors.shape[0] > 0: + FLAs = arrayMeteors[:,4] + maxFLA = FLAs.argmax() + listMeteors1.append(arrayMeteors[maxFLA,:]) - if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'): - theta_x = numpy.array(kwargs['dirCosx']) - theta_y = numpy.array(kwargs['dirCosy']) - else: - elev = numpy.array(kwargs['elevation']) - azim = numpy.array(kwargs['azimuth']) - theta_x, theta_y = self.__calculateCosDir(elev, azim) - azimuth = kwargs['correctAzimuth'] - if kwargs.has_key('horizontalOnly'): - horizontalOnly = kwargs['horizontalOnly'] - else: horizontalOnly = False - if kwargs.has_key('correctFactor'): - correctFactor = kwargs['correctFactor'] - else: correctFactor = 1 - if kwargs.has_key('channelList'): - channelList = kwargs['channelList'] - if len(channelList) == 2: - horizontalOnly = True - arrayChannel = numpy.array(channelList) - param = param[arrayChannel,:,:] - theta_x = theta_x[arrayChannel] - theta_y = theta_y[arrayChannel] + MeteorInitTime = arrayMeteors[maxFLA,1] + MeteorEndTime = arrayMeteors[maxFLA,3] + MeteorHeight = arrayMeteors[maxFLA,0] - velRadial0 = param[:,1,:] #Radial velocity - dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function - dataOut.utctimeInit = dataOut.utctime - dataOut.outputInterval = dataOut.paramInterval + #Check neighborhood + maxHeightIndex = MeteorHeight + rangeLimit + minHeightIndex = MeteorHeight - rangeLimit + minTimeIndex = MeteorInitTime - timeLimit + maxTimeIndex = MeteorEndTime + timeLimit - elif technique == 'SA': + #Check Heights + indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex) + indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex) + indBoth = numpy.where(numpy.logical_and(indTime,indHeight)) + + arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0) - #Parameters - position_x = kwargs['positionX'] - position_y = kwargs['positionY'] - azimuth = kwargs['azimuth'] + return listMeteors1 + + def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency): + numHeights = volts.shape[2] + nChannel = volts.shape[0] + + thresholdPhase = thresh[0] + thresholdNoise = thresh[1] + thresholdDB = float(thresh[2]) + + thresholdDB1 = 10**(thresholdDB/10) + pairsarray = numpy.array(pairslist) + indSides = pairsarray[:,1] + + pairslist1 = list(pairslist) + pairslist1.append((0,1)) + pairslist1.append((3,4)) + + listMeteors1 = [] + listPowerSeries = [] + listVoltageSeries = [] + #volts has the war data + + if frequency == 30e6: + timeLag = 45*10**-3 + else: + timeLag = 15*10**-3 + lag = numpy.ceil(timeLag/timeInterval) + + for i in range(len(listMeteors)): + + ###################### 3.6 - 3.7 PARAMETERS REESTIMATION ######################### + meteorAux = numpy.zeros(16) + + #Loading meteor Data (mHeight, mStart, mPeak, mEnd) + mHeight = listMeteors[i][0] + mStart = listMeteors[i][1] + mPeak = listMeteors[i][2] + mEnd = listMeteors[i][3] + + #get the volt data between the start and end times of the meteor + meteorVolts = volts[:,mStart:mEnd+1,mHeight] + meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) - if kwargs.has_key('crosspairsList'): - pairs = kwargs['crosspairsList'] - else: - pairs = None - - if kwargs.has_key('correctFactor'): - correctFactor = kwargs['correctFactor'] - else: - correctFactor = 1 - - tau = dataOut.data_param - _lambda = dataOut.C/dataOut.frequency - pairsList = dataOut.groupList - nChannels = dataOut.nChannels - - dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) - dataOut.utctimeInit = dataOut.utctime - dataOut.outputInterval = dataOut.timeInterval + #3.6. Phase Difference estimation + phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist) - elif technique == 'Meteors': - dataOut.flagNoData = True - self.__dataReady = False + #3.7. Phase difference removal & meteor start, peak and end times reestimated + #meteorVolts0.- all Channels, all Profiles + meteorVolts0 = volts[:,:,mHeight] + meteorThresh = noise[:,mHeight]*thresholdNoise + meteorNoise = noise[:,mHeight] + meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting + powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power - if kwargs.has_key('nHours'): - nHours = kwargs['nHours'] - else: - nHours = 1 - - if kwargs.has_key('meteorsPerBin'): - meteorThresh = kwargs['meteorsPerBin'] - else: - meteorThresh = 6 + #Times reestimation + mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0] + if mStart1.size > 0: + mStart1 = mStart1[-1] + 1 - if kwargs.has_key('hmin'): - hmin = kwargs['hmin'] - else: hmin = 70 - if kwargs.has_key('hmax'): - hmax = kwargs['hmax'] - else: hmax = 110 - - dataOut.outputInterval = nHours*3600 - - if self.__isConfig == False: -# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) - #Get Initial LTC time - self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) - self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() - - self.__isConfig = True + else: + mStart1 = mPeak - if self.__buffer == None: - self.__buffer = dataOut.data_param - self.__firstdata = copy.copy(dataOut) - + mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1 + mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0] + if mEndDecayTime1.size == 0: + mEndDecayTime1 = powerNet0.size else: - self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1 +# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax() - self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready + #meteorVolts1.- all Channels, from start to end + meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1] + meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1] + if meteorVolts2.shape[1] == 0: + meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1] + meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1) + meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1) + ##################### END PARAMETERS REESTIMATION ######################### - if self.__dataReady: - dataOut.utctimeInit = self.__initime + ##################### 3.8 PHASE DIFFERENCE REESTIMATION ######################## +# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis + if meteorVolts2.shape[1] > 0: + #Phase Difference re-estimation + phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation +# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist) + meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1]) + phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1)) + meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting - self.__initime += dataOut.outputInterval #to erase time offset + #Phase Difference RMS + phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1))) + powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0) + #Data from Meteor + mPeak1 = powerNet1.argmax() + mStart1 + mPeakPower1 = powerNet1.max() + noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight]) + mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux + Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]) + Meteor1 = numpy.hstack((Meteor1,phaseDiffint)) + PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1] + #Vectorize + meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1] + meteorAux[7:11] = phaseDiffint[0:4] - dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) - dataOut.flagNoData = False - self.__buffer = None + #Rejection Criterions + if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation + meteorAux[-1] = 17 + elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB + meteorAux[-1] = 1 - elif technique == 'Meteors1': - dataOut.flagNoData = True - self.__dataReady = False - - if kwargs.has_key('nMins'): - nMins = kwargs['nMins'] - else: nMins = 20 - if kwargs.has_key('rx_location'): - rx_location = kwargs['rx_location'] - else: rx_location = [(0,1),(1,1),(1,0)] - if kwargs.has_key('azimuth'): - azimuth = kwargs['azimuth'] - else: azimuth = 51 - if kwargs.has_key('dfactor'): - dfactor = kwargs['dfactor'] - if kwargs.has_key('mode'): - mode = kwargs['mode'] - else: mode = 'SA' - #Borrar luego esto - if dataOut.groupList == None: - dataOut.groupList = [(0,1),(0,2),(1,2)] - groupList = dataOut.groupList - C = 3e8 - freq = 50e6 - lamb = C/freq - k = 2*numpy.pi/lamb + else: + meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd] + meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis + PowerSeries = 0 + + listMeteors1.append(meteorAux) + listPowerSeries.append(PowerSeries) + listVoltageSeries.append(meteorVolts1) + + return listMeteors1, listPowerSeries, listVoltageSeries + + def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency): + + threshError = 10 + #Depending if it is 30 or 50 MHz + if frequency == 30e6: + timeLag = 45*10**-3 + else: + timeLag = 15*10**-3 + lag = numpy.ceil(timeLag/timeInterval) + + listMeteors1 = [] + + for i in range(len(listMeteors)): + meteorPower = listPower[i] + meteorAux = listMeteors[i] - timeList = dataOut.abscissaList - heightList = dataOut.heightList + if meteorAux[-1] == 0: + + try: + indmax = meteorPower.argmax() + indlag = indmax + lag + + y = meteorPower[indlag:] + x = numpy.arange(0, y.size)*timeLag + + #first guess + a = y[0] + tau = timeLag + #exponential fit + popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau]) + y1 = self.__exponential_function(x, *popt) + #error estimation + error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size)) + + decayTime = popt[1] + riseTime = indmax*timeInterval + meteorAux[11:13] = [decayTime, error] + + #Table items 7, 8 and 11 + if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s + meteorAux[-1] = 7 + elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time + meteorAux[-1] = 8 + if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time + meteorAux[-1] = 11 + + + except: + meteorAux[-1] = 11 + - if self.__isConfig == False: - dataOut.outputInterval = nMins*60 -# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) - #Get Initial LTC time - initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) - minuteAux = initime.minute - minuteNew = int(numpy.floor(minuteAux/nMins)*nMins) - self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + listMeteors1.append(meteorAux) + + return listMeteors1 - self.__isConfig = True + #Exponential Function + + def __exponential_function(self, x, a, tau): + y = a*numpy.exp(-x/tau) + return y + + def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval): + + pairslist1 = list(pairslist) + pairslist1.append((0,1)) + pairslist1.append((3,4)) + numPairs = len(pairslist1) + #Time Lag + timeLag = 45*10**-3 + c = 3e8 + lag = numpy.ceil(timeLag/timeInterval) + freq = 30e6 + + listMeteors1 = [] + + for i in range(len(listMeteors)): + meteorAux = listMeteors[i] + if meteorAux[-1] == 0: + mStart = listMeteors[i][1] + mPeak = listMeteors[i][2] + mLag = mPeak - mStart + lag + + #get the volt data between the start and end times of the meteor + meteorVolts = listVolts[i] + meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) + + #Get CCF + allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2]) + + #Method 2 + slopes = numpy.zeros(numPairs) + time = numpy.array([-2,-1,1,2])*timeInterval + angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0]) + + #Correct phases + derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1] + indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) - if self.__buffer == None: - self.__buffer = dataOut.data_param - self.__firstdata = copy.copy(dataOut) + if indDer[0].shape[0] > 0: + for i in range(indDer[0].shape[0]): + signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]]) + angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi - else: - self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - - self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready - - if self.__dataReady: - dataOut.utctimeInit = self.__initime - self.__initime += dataOut.outputInterval #to erase time offset +# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]])) + for j in range(numPairs): + fit = stats.linregress(time, angAllCCF[j,:]) + slopes[j] = fit[0] - metArray = self.__buffer - if mode == 'SA': - dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList) - elif mode == 'DBS': - dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList) - dataOut.data_output = dataOut.data_output.T - dataOut.flagNoData = False - self.__buffer = None - - return - -class EWDriftsEstimation(Operation): - - def __init__(self): - Operation.__init__(self) + #Remove Outlier +# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) +# slopes = numpy.delete(slopes,indOut) +# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) +# slopes = numpy.delete(slopes,indOut) + + radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq) + radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq) + meteorAux[-2] = radialError + meteorAux[-3] = radialVelocity + + #Setting Error + #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s + if numpy.abs(radialVelocity) > 200: + meteorAux[-1] = 15 + #Number 12: Poor fit to CCF variation for estimation of radial drift velocity + elif radialError > radialStdThresh: + meteorAux[-1] = 12 + + listMeteors1.append(meteorAux) + return listMeteors1 - def __correctValues(self, heiRang, phi, velRadial, SNR): - listPhi = phi.tolist() - maxid = listPhi.index(max(listPhi)) - minid = listPhi.index(min(listPhi)) + def __setNewArrays(self, listMeteors, date, heiRang): - rango = range(len(phi)) - # rango = numpy.delete(rango,maxid) + #New arrays + arrayMeteors = numpy.array(listMeteors) + arrayParameters = numpy.zeros((len(listMeteors), 13)) - heiRang1 = heiRang*math.cos(phi[maxid]) - heiRangAux = heiRang*math.cos(phi[minid]) - indOut = (heiRang1 < heiRangAux[0]).nonzero() - heiRang1 = numpy.delete(heiRang1,indOut) + #Date inclusion +# date = re.findall(r'\((.*?)\)', date) +# date = date[0].split(',') +# date = map(int, date) +# +# if len(date)<6: +# date.append(0) +# +# date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] +# arrayDate = numpy.tile(date, (len(listMeteors), 1)) + arrayDate = numpy.tile(date, (len(listMeteors))) - velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) - SNR1 = numpy.zeros([len(phi),len(heiRang1)]) + #Meteor array +# arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] +# arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) - for i in rango: - x = heiRang*math.cos(phi[i]) - y1 = velRadial[i,:] - f1 = interpolate.interp1d(x,y1,kind = 'cubic') - - x1 = heiRang1 - y11 = f1(x1) - - y2 = SNR[i,:] - f2 = interpolate.interp1d(x,y2,kind = 'cubic') - y21 = f2(x1) - - velRadial1[i,:] = y11 - SNR1[i,:] = y21 - - return heiRang1, velRadial1, SNR1 + #Parameters Array + arrayParameters[:,0] = arrayDate #Date + arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range + arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error + arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases + arrayParameters[:,-1] = arrayMeteors[:,-1] #Error - def run(self, dataOut, zenith, zenithCorrection): - heiRang = dataOut.heightList - velRadial = dataOut.data_param[:,3,:] - SNR = dataOut.data_SNR - zenith = numpy.array(zenith) - zenith -= zenithCorrection - zenith *= numpy.pi/180 - - heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) - - alp = zenith[0] - bet = zenith[1] - - w_w = velRadial1[0,:] - w_e = velRadial1[1,:] - - w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) - u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) + return arrayParameters + +class CorrectMeteorPhases(Operation): + + def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None): + + arrayParameters = dataOut.data_param + pairsList = [] + pairx = (0,1) + pairy = (2,3) + pairsList.append(pairx) + pairsList.append(pairy) + jph = numpy.zeros(4) - winds = numpy.vstack((u,w)) + phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180 + # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets) + arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets))) - dataOut.heightList = heiRang1 - dataOut.data_output = winds - dataOut.data_SNR = SNR1 + meteorOps = MeteorOperations() + if channelPositions == None: + # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T + channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella + + pairslist0, distances = meteorOps.getPhasePairs(channelPositions) + h = (hmin,hmax) + + arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph) - dataOut.utctimeInit = dataOut.utctime - dataOut.outputInterval = dataOut.timeInterval + dataOut.data_param = arrayParameters return - + class PhaseCalibration(Operation): __buffer = None @@ -2709,4 +2615,108 @@ class MeteorOperations(): # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)] pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)] - return pairslist, distances \ No newline at end of file + return pairslist, distances +# def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): +# +# arrayAOA = numpy.zeros((phases.shape[0],3)) +# cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList) +# +# arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth) +# cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1) +# arrayAOA[:,2] = cosDirError +# +# azimuthAngle = arrayAOA[:,0] +# zenithAngle = arrayAOA[:,1] +# +# #Setting Error +# #Number 3: AOA not fesible +# indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0] +# error[indInvalid] = 3 +# #Number 4: Large difference in AOAs obtained from different antenna baselines +# indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0] +# error[indInvalid] = 4 +# return arrayAOA, error +# +# def __getDirectionCosines(self, arrayPhase, pairsList): +# +# #Initializing some variables +# ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi +# ang_aux = ang_aux.reshape(1,ang_aux.size) +# +# cosdir = numpy.zeros((arrayPhase.shape[0],2)) +# cosdir0 = numpy.zeros((arrayPhase.shape[0],2)) +# +# +# for i in range(2): +# #First Estimation +# phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]] +# #Dealias +# indcsi = numpy.where(phi0_aux > numpy.pi) +# phi0_aux[indcsi] -= 2*numpy.pi +# indcsi = numpy.where(phi0_aux < -numpy.pi) +# phi0_aux[indcsi] += 2*numpy.pi +# #Direction Cosine 0 +# cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5) +# +# #Most-Accurate Second Estimation +# phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]] +# phi1_aux = phi1_aux.reshape(phi1_aux.size,1) +# #Direction Cosine 1 +# cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5) +# +# #Searching the correct Direction Cosine +# cosdir0_aux = cosdir0[:,i] +# cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) +# #Minimum Distance +# cosDiff = (cosdir1 - cosdir0_aux)**2 +# indcos = cosDiff.argmin(axis = 1) +# #Saving Value obtained +# cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos] +# +# return cosdir0, cosdir +# +# def __calculateAOA(self, cosdir, azimuth): +# cosdirX = cosdir[:,0] +# cosdirY = cosdir[:,1] +# +# zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi +# azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east +# angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose() +# +# return angles +# +# def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight): +# +# Ramb = 375 #Ramb = c/(2*PRF) +# Re = 6371 #Earth Radius +# heights = numpy.zeros(Ranges.shape) +# +# R_aux = numpy.array([0,1,2])*Ramb +# R_aux = R_aux.reshape(1,R_aux.size) +# +# Ranges = Ranges.reshape(Ranges.size,1) +# +# Ri = Ranges + R_aux +# hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re +# +# #Check if there is a height between 70 and 110 km +# h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1) +# ind_h = numpy.where(h_bool == 1)[0] +# +# hCorr = hi[ind_h, :] +# ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight)) +# +# hCorr = hi[ind_hCorr] +# heights[ind_h] = hCorr +# +# #Setting Error +# #Number 13: Height unresolvable echo: not valid height within 70 to 110 km +# #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km +# +# indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] +# error[indInvalid2] = 14 +# indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0] +# error[indInvalid1] = 13 +# +# return heights, error + \ No newline at end of file