##// END OF EJS Templates
schain-jc
RSS
Bug Fixed: AMISR Setup File does not correspond with the experiment range dates...
Bug Fixed: AMISR Setup File does not correspond with the experiment range dates Bug Fixed: Coherent Integration, the input parameter 'timeInterval' should be in seconds unit
Julio Valdez - - Load All Authors

File last commit:

r503:ea980322a40d
r510:de96659ace58
Show More
jroproc_parameters.py
1538 lines | 62.4 KiB | text/x-python | PythonLexer
/ schainpy / model / proc / jroproc_parameters.py
History | Annotation | Raw |Copy content |Copy permalink
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
from jroproc_base import ProcessingUnit, Operation
from model.data.jrodata import Parameters
class ParametersProc(ProcessingUnit):
nSeconds = None
def __init__(self):
ProcessingUnit.__init__(self)
self.objectDict = {}
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
self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
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
self.dataOut.heightRange = self.dataIn.getHeiRange()
self.dataOut.frequency = self.dataIn.frequency
def run(self, nSeconds = None, nProfiles = None):
self.dataOut.flagNoData = True
if self.firstdatatime == None:
self.firstdatatime = self.dataIn.utctime
#---------------------- Voltage Data ---------------------------
if self.dataIn.type == "Voltage":
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
#---------------------- Spectra Data ---------------------------
if self.dataIn.type == "Spectra":
self.dataOut.data_pre = self.dataIn.data_spc.copy()
self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
self.dataOut.noise = self.dataIn.getNoise()
self.dataOut.normFactor = self.dataIn.normFactor
self.__updateObjFromInput()
self.dataOut.flagNoData = False
self.firstdatatime = None
#---------------------- 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,:]
self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
self.dataOut.noise = self.dataIn.noise
self.dataOut.normFactor = self.dataIn.normFactor
self.dataOut.SNR = self.dataIn.SNR
self.dataOut.pairsList = self.dataIn.pairsList
self.__updateObjFromInput()
self.dataOut.flagNoData = False
self.firstdatatime = None
#------------------- 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
self.dataOut.abscissaRange
self.dataOut.noise
Affected:
self.dataOut.data_param
self.dataOut.SNR
'''
data = self.dataOut.data_pre
absc = self.dataOut.abscissaRange[:-1]
noise = self.dataOut.noise
data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
if channelList== None: channelList = self.dataOut.channelList
for ind in channelList:
data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
self.dataOut.data_param = data_param[:,1:]
self.dataOut.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):
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
#------------------- Get Lags ----------------------------------
def GetLags(self):
'''
Function GetMoments()
Input:
self.dataOut.data_pre
self.dataOut.abscissaRange
self.dataOut.noise
self.dataOut.normFactor
self.dataOut.SNR
self.dataOut.pairsList
self.dataOut.nChannels
Affected:
self.dataOut.data_param
'''
data = self.dataOut.data_pre
normFactor = self.dataOut.normFactor
nHeights = self.dataOut.nHeights
absc = self.dataOut.abscissaRange[:-1]
noise = self.dataOut.noise
SNR = self.dataOut.SNR
pairsList = self.dataOut.pairsList
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
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]]
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
def __checkTime(self, currentTime, paramInterval, windsInterval):
dataTime = currentTime + paramInterval
deltaTime = dataTime - self.__initime
if deltaTime >= windsInterval or deltaTime < 0:
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
if dataOut.abscissaRange != None:
absc = dataOut.abscissaRange[:-1]
noise = dataOut.noise
heightRange = dataOut.getHeiRange()
SNR = dataOut.SNR
if technique == 'DBS':
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]
velRadial0 = param[:,1,:] #Radial velocity
dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
dataOut.initUtcTime = dataOut.ltctime
dataOut.windsInterval = dataOut.timeInterval
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
pairsList = dataOut.pairsList
nChannels = dataOut.nChannels
dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
dataOut.initUtcTime = dataOut.ltctime
dataOut.windsInterval = dataOut.timeInterval
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
dataOut.windsInterval = nHours*3600
if self.__isConfig == False:
# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
#Get Initial LTC time
self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
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))
self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
if self.__dataReady:
dataOut.initUtcTime = self.__initime
self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
dataOut.flagNoData = False
self.__buffer = None
return