jroproc_parameters.py
1779 lines
| 70.7 KiB
| text/x-python
|
PythonLexer
Julio Valdez
|
r502 | import numpy | |
| import math | |||
| from scipy import optimize | |||
| from scipy import interpolate | |||
| from scipy import signal | |||
| from scipy import stats | |||
| import re | |||
| import datetime | |||
| import copy | |||
|
Julio Valdez
|
r513 | import sys | |
| import importlib | |||
| import itertools | |||
|
Julio Valdez
|
r502 | ||
| from jroproc_base import ProcessingUnit, Operation | |||
Miguel Valdez
|
r568 | from schainpy.model.data.jrodata import Parameters | |
|
Julio Valdez
|
r502 | ||
| class ParametersProc(ProcessingUnit): | |||
| nSeconds = None | |||
| def __init__(self): | |||
| ProcessingUnit.__init__(self) | |||
|
Julio Valdez
|
r543 | # self.objectDict = {} | |
|
Julio Valdez
|
r502 | self.buffer = None | |
| self.firstdatatime = None | |||
| self.profIndex = 0 | |||
| self.dataOut = Parameters() | |||
| def __updateObjFromInput(self): | |||
| self.dataOut.inputUnit = self.dataIn.type | |||
| self.dataOut.timeZone = self.dataIn.timeZone | |||
| self.dataOut.dstFlag = self.dataIn.dstFlag | |||
| self.dataOut.errorCount = self.dataIn.errorCount | |||
| self.dataOut.useLocalTime = self.dataIn.useLocalTime | |||
| self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() | |||
| self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy() | |||
| self.dataOut.channelList = self.dataIn.channelList | |||
| self.dataOut.heightList = self.dataIn.heightList | |||
| self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')]) | |||
| # self.dataOut.nHeights = self.dataIn.nHeights | |||
| # self.dataOut.nChannels = self.dataIn.nChannels | |||
| self.dataOut.nBaud = self.dataIn.nBaud | |||
| self.dataOut.nCode = self.dataIn.nCode | |||
| self.dataOut.code = self.dataIn.code | |||
| # self.dataOut.nProfiles = self.dataOut.nFFTPoints | |||
|
Miguel Valdez
|
r568 | self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock | |
|
Julio Valdez
|
r502 | self.dataOut.utctime = self.firstdatatime | |
| self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada | |||
| self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip | |||
| # self.dataOut.nCohInt = self.dataIn.nCohInt | |||
| # self.dataOut.nIncohInt = 1 | |||
| self.dataOut.ippSeconds = self.dataIn.ippSeconds | |||
| # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter | |||
| self.dataOut.timeInterval = self.dataIn.timeInterval | |||
|
Julio Valdez
|
r543 | self.dataOut.heightList = self.dataIn.getHeiRange() | |
|
Julio Valdez
|
r502 | self.dataOut.frequency = self.dataIn.frequency | |
| def run(self, nSeconds = None, nProfiles = None): | |||
|
Julio Valdez
|
r513 | ||
|
Julio Valdez
|
r502 | ||
| if self.firstdatatime == None: | |||
| self.firstdatatime = self.dataIn.utctime | |||
| #---------------------- Voltage Data --------------------------- | |||
| if self.dataIn.type == "Voltage": | |||
|
Julio Valdez
|
r513 | self.dataOut.flagNoData = True | |
|
Julio Valdez
|
r502 | if nSeconds != None: | |
| self.nSeconds = nSeconds | |||
| self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt))) | |||
| if self.buffer == None: | |||
| self.buffer = numpy.zeros((self.dataIn.nChannels, | |||
| self.nProfiles, | |||
| self.dataIn.nHeights), | |||
| dtype='complex') | |||
| self.buffer[:,self.profIndex,:] = self.dataIn.data.copy() | |||
| self.profIndex += 1 | |||
| if self.profIndex == self.nProfiles: | |||
| self.__updateObjFromInput() | |||
| self.dataOut.data_pre = self.buffer.copy() | |||
| self.dataOut.paramInterval = nSeconds | |||
| self.dataOut.flagNoData = False | |||
| self.buffer = None | |||
| self.firstdatatime = None | |||
| self.profIndex = 0 | |||
|
Julio Valdez
|
r511 | return | |
|
Julio Valdez
|
r502 | ||
| #---------------------- Spectra Data --------------------------- | |||
| if self.dataIn.type == "Spectra": | |||
| self.dataOut.data_pre = self.dataIn.data_spc.copy() | |||
|
Julio Valdez
|
r543 | self.dataOut.abscissaList = self.dataIn.getVelRange(1) | |
|
Julio Valdez
|
r502 | self.dataOut.noise = self.dataIn.getNoise() | |
| self.dataOut.normFactor = self.dataIn.normFactor | |||
|
Julio Valdez
|
r549 | self.dataOut.groupList = self.dataIn.pairsList | |
|
Julio Valdez
|
r513 | self.dataOut.flagNoData = False | |
|
Julio Valdez
|
r511 | ||
|
Julio Valdez
|
r502 | #---------------------- Correlation Data --------------------------- | |
| if self.dataIn.type == "Correlation": | |||
| lagRRange = self.dataIn.lagR | |||
| indR = numpy.where(lagRRange == 0)[0][0] | |||
| self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:] | |||
|
Julio Valdez
|
r543 | self.dataOut.abscissaList = self.dataIn.getLagTRange(1) | |
|
Julio Valdez
|
r502 | self.dataOut.noise = self.dataIn.noise | |
| self.dataOut.normFactor = self.dataIn.normFactor | |||
|
Julio Valdez
|
r514 | self.dataOut.data_SNR = self.dataIn.SNR | |
|
Julio Valdez
|
r513 | self.dataOut.groupList = self.dataIn.pairsList | |
|
Julio Valdez
|
r543 | self.dataOut.flagNoData = False | |
| #---------------------- Correlation Data --------------------------- | |||
| if self.dataIn.type == "Parameters": | |||
| self.dataOut.copy(self.dataIn) | |||
|
Julio Valdez
|
r513 | self.dataOut.flagNoData = False | |
|
Julio Valdez
|
r502 | ||
|
Julio Valdez
|
r543 | return True | |
|
Julio Valdez
|
r511 | ||
| self.__updateObjFromInput() | |||
| self.firstdatatime = None | |||
|
Julio Valdez
|
r543 | self.dataOut.utctimeInit = self.dataIn.utctime | |
|
Julio Valdez
|
r513 | self.dataOut.outputInterval = self.dataIn.timeInterval | |
|
Julio Valdez
|
r502 | ||
| #------------------- Get Moments ---------------------------------- | |||
| def GetMoments(self, channelList = None): | |||
| ''' | |||
| Function GetMoments() | |||
| Input: | |||
| channelList : simple channel list to select e.g. [2,3,7] | |||
| self.dataOut.data_pre | |||
|
Julio Valdez
|
r543 | self.dataOut.abscissaList | |
|
Julio Valdez
|
r502 | self.dataOut.noise | |
| Affected: | |||
| self.dataOut.data_param | |||
|
Julio Valdez
|
r514 | self.dataOut.data_SNR | |
|
Julio Valdez
|
r502 | ||
| ''' | |||
| data = self.dataOut.data_pre | |||
|
Julio Valdez
|
r543 | absc = self.dataOut.abscissaList[:-1] | |
|
Julio Valdez
|
r502 | noise = self.dataOut.noise | |
| data_param = numpy.zeros((data.shape[0], 4, data.shape[2])) | |||
|
Julio Valdez
|
r511 | if channelList== None: | |
| channelList = self.dataIn.channelList | |||
| self.dataOut.channelList = channelList | |||
|
Julio Valdez
|
r502 | ||
| for ind in channelList: | |||
| data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind]) | |||
|
Julio Valdez
|
r511 | ||
| self.dataOut.data_param = data_param[:,1:,:] | |||
|
Julio Valdez
|
r514 | self.dataOut.data_SNR = data_param[:,0] | |
|
Julio Valdez
|
r502 | 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): | |||
| if (nicoh == None): nicoh = 1 | |||
| if (graph == None): graph = 0 | |||
| if (smooth == None): smooth = 0 | |||
| elif (self.smooth < 3): smooth = 0 | |||
| if (type1 == None): type1 = 0 | |||
| if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1 | |||
| if (snrth == None): snrth = -3 | |||
| if (dc == None): dc = 0 | |||
| if (aliasing == None): aliasing = 0 | |||
| if (oldfd == None): oldfd = 0 | |||
| if (wwauto == None): wwauto = 0 | |||
| if (n0 < 1.e-20): n0 = 1.e-20 | |||
| freq = oldfreq | |||
| vec_power = numpy.zeros(oldspec.shape[1]) | |||
| vec_fd = numpy.zeros(oldspec.shape[1]) | |||
| vec_w = numpy.zeros(oldspec.shape[1]) | |||
| vec_snr = numpy.zeros(oldspec.shape[1]) | |||
| for ind in range(oldspec.shape[1]): | |||
| spec = oldspec[:,ind] | |||
| aux = spec*fwindow | |||
| max_spec = aux.max() | |||
| m = list(aux).index(max_spec) | |||
| #Smooth | |||
| if (smooth == 0): spec2 = spec | |||
| else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth) | |||
| # Calculo de Momentos | |||
| bb = spec2[range(m,spec2.size)] | |||
| bb = (bb<n0).nonzero() | |||
| bb = bb[0] | |||
| ss = spec2[range(0,m + 1)] | |||
| ss = (ss<n0).nonzero() | |||
| ss = ss[0] | |||
| if (bb.size == 0): | |||
| bb0 = spec.size - 1 - m | |||
| else: | |||
| bb0 = bb[0] - 1 | |||
| if (bb0 < 0): | |||
| bb0 = 0 | |||
| if (ss.size == 0): ss1 = 1 | |||
| else: ss1 = max(ss) + 1 | |||
| if (ss1 > m): ss1 = m | |||
| valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1 | |||
| power = ((spec2[valid] - n0)*fwindow[valid]).sum() | |||
| fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power | |||
| w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) | |||
| snr = (spec2.mean()-n0)/n0 | |||
| if (snr < 1.e-20) : | |||
| snr = 1.e-20 | |||
| vec_power[ind] = power | |||
| vec_fd[ind] = fd | |||
| vec_w[ind] = w | |||
| vec_snr[ind] = snr | |||
| moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) | |||
| return moments | |||
|
Julio Valdez
|
r549 | #------------------ Get SA Parameters -------------------------- | |
| def GetSAParameters(self): | |||
| data = self.dataOut.data_pre | |||
| crossdata = self.dataIn.data_cspc | |||
| a = 1 | |||
|
Julio Valdez
|
r502 | #------------------- Get Lags ---------------------------------- | |
| def GetLags(self): | |||
| ''' | |||
| Function GetMoments() | |||
| Input: | |||
| self.dataOut.data_pre | |||
|
Julio Valdez
|
r543 | self.dataOut.abscissaList | |
|
Julio Valdez
|
r502 | self.dataOut.noise | |
| self.dataOut.normFactor | |||
|
Julio Valdez
|
r514 | self.dataOut.data_SNR | |
|
Julio Valdez
|
r513 | self.dataOut.groupList | |
|
Julio Valdez
|
r502 | self.dataOut.nChannels | |
| Affected: | |||
| self.dataOut.data_param | |||
| ''' | |||
|
Julio Valdez
|
r549 | ||
|
Julio Valdez
|
r502 | data = self.dataOut.data_pre | |
| normFactor = self.dataOut.normFactor | |||
| nHeights = self.dataOut.nHeights | |||
|
Julio Valdez
|
r543 | absc = self.dataOut.abscissaList[:-1] | |
|
Julio Valdez
|
r502 | noise = self.dataOut.noise | |
|
Julio Valdez
|
r514 | SNR = self.dataOut.data_SNR | |
|
Julio Valdez
|
r513 | pairsList = self.dataOut.groupList | |
|
Julio Valdez
|
r502 | nChannels = self.dataOut.nChannels | |
| pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) | |||
| self.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) | |||
| 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 __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange): | |||
| Pt0 = data.shape[1]/2 | |||
| #Funcion de Autocorrelacion | |||
| dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0) | |||
| #Obtencion Indice de TauCross | |||
| indCross = data[pairsCrossCorr,:,:].argmax(axis = 1) | |||
| #Obtencion Indice de TauAuto | |||
| indAuto = numpy.zeros(indCross.shape,dtype = 'int') | |||
| CCValue = data[pairsCrossCorr,Pt0,:] | |||
| for i in range(pairsCrossCorr.size): | |||
| indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0) | |||
| #Obtencion de TauCross y TauAuto | |||
| tauCross = lagTRange[indCross] | |||
| tauAuto = lagTRange[indAuto] | |||
| Nan1, Nan2 = numpy.where(tauCross == lagTRange[0]) | |||
| tauCross[Nan1,Nan2] = numpy.nan | |||
| tauAuto[Nan1,Nan2] = numpy.nan | |||
| tau = numpy.vstack((tauCross,tauAuto)) | |||
| return tau | |||
| def __calculateLag1Phase(self, data, pairs, lagTRange): | |||
| data1 = stats.nanmean(data[pairs,:,:], axis = 0) | |||
| lag1 = numpy.where(lagTRange == 0)[0][0] + 1 | |||
| phase = numpy.angle(data1[lag1,:]) | |||
| return phase | |||
| #------------------- Detect Meteors ------------------------------ | |||
| def DetectMeteors(self, hei_ref = None, tauindex = 0, | |||
| predefinedPhaseShifts = None, centerReceiverIndex = 2, | |||
| cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, | |||
| noise_timeStep = 4, noise_multiple = 4, | |||
| multDet_timeLimit = 1, multDet_rangeLimit = 3, | |||
| phaseThresh = 20, SNRThresh = 8, | |||
| hmin = 70, hmax=110, azimuth = 0) : | |||
| ''' | |||
| 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): | |||
| Day Hour | Range Height | |||
| Azimuth Zenith errorCosDir | |||
| VelRad errorVelRad | |||
| TypeError | |||
| ''' | |||
| #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) | |||
| heiRang = self.dataOut.getHeiRange() | |||
| #Pairs List | |||
| pairslist = [] | |||
| 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() | |||
| if predefinedPhaseShifts != None: | |||
| hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 | |||
| else: | |||
| #get hardware phase shifts using beacon signal | |||
| hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) | |||
| hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) | |||
| 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]] | |||
| #************ 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, pairslist, cohDet_thresh) | |||
| #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 ********** | |||
| #************** 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 ********************** | |||
| #********************* 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, pairslist, 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, pairslist, self.dataOut.timeInterval) | |||
| if len(listMeteors4) > 0: | |||
| #Setting New Array | |||
| date = repr(self.dataOut.datatime) | |||
| arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) | |||
| #Calculate AOA (Error N 3, 4) | |||
| #JONES ET AL. 1998 | |||
| AOAthresh = numpy.pi/8 | |||
| error = arrayParameters[:,-1] | |||
| phases = -arrayMeteors4[:,9:13] | |||
| pairsList = [] | |||
| pairsList.append((0,3)) | |||
| pairsList.append((1,2)) | |||
| arrayParameters[:,4:7], arrayParameters[:,-1] = self.__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] = self.__getHeights(Ranges, zenith, error, hmin, hmax) | |||
| #********************* END OF PARAMETERS CALCULATION ************************** | |||
| #***************************+ SAVE DATA IN HDF5 FORMAT ********************** | |||
| self.dataOut.data_param = arrayParameters | |||
| return | |||
| 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) | |||
| return phaseOffset | |||
| 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])) | |||
| #Correct phases | |||
| derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:] | |||
| indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) | |||
| 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 | |||
| # 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) | |||
| #Dealias | |||
| indAlias = numpy.where(phaseArrival > numpy.pi) | |||
| phaseArrival[indAlias] -= 2*numpy.pi | |||
| indAlias = numpy.where(phaseArrival < -numpy.pi) | |||
| phaseArrival[indAlias] += 2*numpy.pi | |||
| 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 | |||
| 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): | |||
| 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) | |||
| vStacked = None | |||
| return voltsCCF | |||
| 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])) | |||
| startInd = 0 | |||
| endInd = 0 | |||
| 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 = [] | |||
| 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 | |||
| def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit): | |||
| arrayMeteors = numpy.asarray(listMeteors) | |||
| listMeteors1 = [] | |||
| 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)) | |||
| arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0) | |||
| 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) | |||
| #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): | |||
| 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] | |||
| 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) | |||
| return listMeteors1 | |||
| #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)): | |||
| meteor = listMeteors[i] | |||
| meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1])) | |||
| if meteor[-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 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] | |||
| #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 __setNewArrays(self, listMeteors, date, heiRang): | |||
| #New arrays | |||
| arrayMeteors = numpy.array(listMeteors) | |||
| arrayParameters = numpy.zeros((len(listMeteors),10)) | |||
| #Date inclusion | |||
| date = re.findall(r'\((.*?)\)', date) | |||
| date = date[0].split(',') | |||
| date = map(int, date) | |||
| date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] | |||
| arrayDate = numpy.tile(date, (len(listMeteors), 1)) | |||
| #Meteor array | |||
| arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] | |||
| arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) | |||
| #Parameters Array | |||
| arrayParameters[:,0:3] = arrayMeteors[:,0:3] | |||
| arrayParameters[:,-3:] = arrayMeteors[:,-3:] | |||
| return arrayMeteors, 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 | |||
|
Julio Valdez
|
r513 | def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None): | |
| ''' | |||
| Function GetMoments() | |||
| Input: | |||
| Output: | |||
| Variables modified: | |||
| ''' | |||
| if path != None: | |||
| sys.path.append(path) | |||
| self.dataOut.library = importlib.import_module(file) | |||
| #To be inserted as a parameter | |||
| groupArray = numpy.array(groupList) | |||
| # groupArray = numpy.array([[0,1],[2,3]]) | |||
| self.dataOut.groupList = groupArray | |||
| nGroups = groupArray.shape[0] | |||
| nChannels = self.dataIn.nChannels | |||
| nHeights=self.dataIn.heightList.size | |||
| #Parameters Array | |||
| self.dataOut.data_param = None | |||
| #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) | |||
| #List of possible combinations | |||
| listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) | |||
| indCross = numpy.zeros(len(list(listComb)), dtype = 'int') | |||
| if getSNR: | |||
| listChannels = groupArray.reshape((groupArray.size)) | |||
| listChannels.sort() | |||
| noise = self.dataIn.getNoise() | |||
|
Julio Valdez
|
r514 | self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) | |
|
Julio Valdez
|
r513 | ||
| for i in range(nGroups): | |||
| coord = groupArray[i,:] | |||
| #Input data array | |||
| data = self.dataIn.data_spc[coord,:,:]/(M*N) | |||
| data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) | |||
| #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 | |||
| for h in range(nHeights): | |||
| # print self.dataOut.heightList[h] | |||
| #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] | |||
|
Julio Valdez
|
r543 | ||
| 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)) | |||
|
Julio Valdez
|
r513 | ||
|
Julio Valdez
|
r514 | 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 | |||
|
Julio Valdez
|
r513 | #Save | |
| if self.dataOut.data_param == None: | |||
|
Julio Valdez
|
r514 | 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 | |||
|
Julio Valdez
|
r513 | ||
|
Julio Valdez
|
r514 | self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1)) | |
|
Julio Valdez
|
r513 | self.dataOut.data_param[i,:,h] = minp | |
| return | |||
| def __residFunction(self, p, dp, LT, constants): | |||
| fm = self.dataOut.library.modelFunction(p, constants) | |||
| fmp=numpy.dot(LT,fm) | |||
| return dp-fmp | |||
| 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 | |||
|
Julio Valdez
|
r502 | ||
| 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) | |||
|
Julio Valdez
|
r503 | 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 | |||
|
Julio Valdez
|
r502 | 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]] | |||
| return distx,disty, dist1,ang1 | |||
| def __calculateVelVer(self, phase, lagTRange, _lambda): | |||
| 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])) | |||
| angCos = numpy.cos(ang) | |||
| angSin = numpy.sin(ang) | |||
| vel0 = dist*tau1/(2*tau2**2) | |||
| vel[:,0,:] = (vel0*angCos).sum(axis = 1) | |||
| vel[:,1,:] = (vel0*angSin).sum(axis = 1) | |||
| 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 | |||
| for l in range(len(pairsList)): | |||
| firstChannel = pairsList[l][0] | |||
| secondChannel = pairsList[l][1] | |||
| #Obteniendo pares de Autocorrelacion | |||
| if firstChannel == secondChannel: | |||
| pairsAutoCorr[firstChannel] = int(l) | |||
| pairsAutoCorr = pairsAutoCorr.astype(int) | |||
| pairsCrossCorr = range(len(pairsList)) | |||
| pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) | |||
| return pairsAutoCorr, pairsCrossCorr | |||
| def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): | |||
| """ | |||
| 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 | |||
| Output: Winds estimation (Zonal, Meridional and Vertical) | |||
| Parameters affected: Winds | |||
| """ | |||
| #Cross Correlation pairs obtained | |||
| pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) | |||
| pairsArray = numpy.array(pairsList)[pairsCrossCorr] | |||
| pairsSelArray = numpy.array(pairsSelected) | |||
| pairs = [] | |||
| #Wind estimation pairs obtained | |||
| for i in range(pairsSelArray.shape[0]/2): | |||
| ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0] | |||
| ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0] | |||
| pairs.append((ind1,ind2)) | |||
| indtau = tau.shape[0]/2 | |||
| tau1 = tau[:indtau,:] | |||
| tau2 = tau[indtau:-1,:] | |||
| tau1 = tau1[pairs,:] | |||
| tau2 = tau2[pairs,:] | |||
| phase1 = tau[-1,:] | |||
| #--------------------------------------------------------------------- | |||
| #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) | |||
| #--------------------------------------------------------------------- | |||
| winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) | |||
| winds = correctFactor*winds | |||
| return winds | |||
|
Julio Valdez
|
r513 | def __checkTime(self, currentTime, paramInterval, outputInterval): | |
|
Julio Valdez
|
r502 | ||
| dataTime = currentTime + paramInterval | |||
| deltaTime = dataTime - self.__initime | |||
|
Julio Valdez
|
r513 | if deltaTime >= outputInterval or deltaTime < 0: | |
|
Julio Valdez
|
r502 | self.__dataReady = True | |
| return | |||
| 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 | |||
| ''' | |||
| #Settings | |||
| nInt = (heightMax - heightMin)/2 | |||
| winds = numpy.zeros((2,nInt))*numpy.nan | |||
| #Filter errors | |||
| error = numpy.where(arrayMeteor[:,-1] == 0)[0] | |||
| finalMeteor = arrayMeteor[error,:] | |||
| #Meteor Histogram | |||
| finalHeights = finalMeteor[:,3] | |||
| 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[:, 7] | |||
| zen = meteorAux[:, 5]*numpy.pi/180 | |||
| azim = meteorAux[:, 4]*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] | |||
| return winds, heightPerI[:-1] | |||
| def run(self, dataOut, technique, **kwargs): | |||
| param = dataOut.data_param | |||
|
Julio Valdez
|
r543 | if dataOut.abscissaList != None: | |
| absc = dataOut.abscissaList[:-1] | |||
|
Julio Valdez
|
r502 | noise = dataOut.noise | |
|
Julio Valdez
|
r543 | heightList = dataOut.getHeiRange() | |
|
Julio Valdez
|
r514 | SNR = dataOut.data_SNR | |
|
Julio Valdez
|
r502 | ||
| if technique == 'DBS': | |||
|
Julio Valdez
|
r503 | 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'] | |||
|
Julio Valdez
|
r502 | 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] | |||
| velRadial0 = param[:,1,:] #Radial velocity | |||
|
Julio Valdez
|
r543 | 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.timeInterval | |||
|
Julio Valdez
|
r502 | ||
| elif technique == 'SA': | |||
| #Parameters | |||
| position_x = kwargs['positionX'] | |||
| position_y = kwargs['positionY'] | |||
| azimuth = kwargs['azimuth'] | |||
| 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 | |||
|
Julio Valdez
|
r513 | pairsList = dataOut.groupList | |
|
Julio Valdez
|
r502 | nChannels = dataOut.nChannels | |
|
Julio Valdez
|
r513 | dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) | |
|
Julio Valdez
|
r543 | dataOut.utctimeInit = dataOut.utctime | |
|
Julio Valdez
|
r513 | dataOut.outputInterval = dataOut.timeInterval | |
|
Julio Valdez
|
r502 | ||
| elif technique == 'Meteors': | |||
| dataOut.flagNoData = True | |||
| self.__dataReady = False | |||
| if kwargs.has_key('nHours'): | |||
| nHours = kwargs['nHours'] | |||
| else: | |||
| nHours = 1 | |||
| if kwargs.has_key('meteorsPerBin'): | |||
| meteorThresh = kwargs['meteorsPerBin'] | |||
| else: | |||
| meteorThresh = 6 | |||
| if kwargs.has_key('hmin'): | |||
| hmin = kwargs['hmin'] | |||
| else: hmin = 70 | |||
| if kwargs.has_key('hmax'): | |||
| hmax = kwargs['hmax'] | |||
| else: hmax = 110 | |||
|
Julio Valdez
|
r513 | dataOut.outputInterval = nHours*3600 | |
|
Julio Valdez
|
r502 | ||
| if self.__isConfig == False: | |||
| # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) | |||
| #Get Initial LTC time | |||
|
Julio Valdez
|
r543 | self.__initime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime) | |
| self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() | |||
|
Julio Valdez
|
r502 | self.__isConfig = True | |
| if self.__buffer == None: | |||
| self.__buffer = dataOut.data_param | |||
| self.__firstdata = copy.copy(dataOut) | |||
| else: | |||
| self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) | |||
|
Julio Valdez
|
r543 | self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready | |
|
Julio Valdez
|
r502 | ||
| if self.__dataReady: | |||
|
Julio Valdez
|
r543 | dataOut.utctimeInit = self.__initime | |
| self.__initime += dataOut.outputInterval #to erase time offset | |||
|
Julio Valdez
|
r502 | ||
|
Julio Valdez
|
r543 | dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) | |
|
Julio Valdez
|
r502 | dataOut.flagNoData = False | |
| self.__buffer = None | |||
|
Julio Valdez
|
r513 | 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,:] | |||
|
Julio Valdez
|
r514 | SNR = dataOut.data_SNR | |
|
Julio Valdez
|
r513 | ||
| 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 | |||
|
Julio Valdez
|
r514 | dataOut.data_SNR = SNR1 | |
|
Julio Valdez
|
r513 | ||
|
Julio Valdez
|
r543 | dataOut.utctimeInit = dataOut.utctime | |
|
Julio Valdez
|
r513 | dataOut.outputInterval = dataOut.timeInterval | |
| return | |||