|
|
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
|
|
|
import sys
|
|
|
import importlib
|
|
|
import itertools
|
|
|
|
|
|
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.heightList = self.dataIn.getHeiRange()
|
|
|
self.dataOut.frequency = self.dataIn.frequency
|
|
|
|
|
|
def run(self, nSeconds = None, nProfiles = None):
|
|
|
|
|
|
|
|
|
|
|
|
if self.firstdatatime == None:
|
|
|
self.firstdatatime = self.dataIn.utctime
|
|
|
|
|
|
#---------------------- Voltage Data ---------------------------
|
|
|
|
|
|
if self.dataIn.type == "Voltage":
|
|
|
self.dataOut.flagNoData = True
|
|
|
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
|
|
|
return
|
|
|
|
|
|
#---------------------- Spectra Data ---------------------------
|
|
|
|
|
|
if self.dataIn.type == "Spectra":
|
|
|
self.dataOut.data_pre = self.dataIn.data_spc.copy()
|
|
|
self.dataOut.abscissaList = self.dataIn.getVelRange(1)
|
|
|
self.dataOut.noise = self.dataIn.getNoise()
|
|
|
self.dataOut.normFactor = self.dataIn.normFactor
|
|
|
self.dataOut.groupList = self.dataIn.pairsList
|
|
|
self.dataOut.flagNoData = False
|
|
|
|
|
|
#---------------------- 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.abscissaList = self.dataIn.getLagTRange(1)
|
|
|
self.dataOut.noise = self.dataIn.noise
|
|
|
self.dataOut.normFactor = self.dataIn.normFactor
|
|
|
self.dataOut.data_SNR = self.dataIn.SNR
|
|
|
self.dataOut.groupList = self.dataIn.pairsList
|
|
|
self.dataOut.flagNoData = False
|
|
|
|
|
|
#---------------------- Correlation Data ---------------------------
|
|
|
|
|
|
if self.dataIn.type == "Parameters":
|
|
|
self.dataOut.copy(self.dataIn)
|
|
|
self.dataOut.flagNoData = False
|
|
|
|
|
|
return True
|
|
|
|
|
|
self.__updateObjFromInput()
|
|
|
self.firstdatatime = None
|
|
|
self.dataOut.utctimeInit = self.dataIn.utctime
|
|
|
self.dataOut.outputInterval = self.dataIn.timeInterval
|
|
|
|
|
|
#------------------- 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.abscissaList
|
|
|
self.dataOut.noise
|
|
|
|
|
|
Affected:
|
|
|
self.dataOut.data_param
|
|
|
self.dataOut.data_SNR
|
|
|
|
|
|
'''
|
|
|
data = self.dataOut.data_pre
|
|
|
absc = self.dataOut.abscissaList[:-1]
|
|
|
noise = self.dataOut.noise
|
|
|
|
|
|
data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
|
|
|
|
|
|
if channelList== None:
|
|
|
channelList = self.dataIn.channelList
|
|
|
self.dataOut.channelList = channelList
|
|
|
|
|
|
for ind in channelList:
|
|
|
data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
|
|
|
|
|
|
self.dataOut.data_param = data_param[:,1:,:]
|
|
|
self.dataOut.data_SNR = data_param[:,0]
|
|
|
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 SA Parameters --------------------------
|
|
|
def GetSAParameters(self):
|
|
|
data = self.dataOut.data_pre
|
|
|
crossdata = self.dataIn.data_cspc
|
|
|
a = 1
|
|
|
|
|
|
|
|
|
|
|
|
#------------------- Get Lags ----------------------------------
|
|
|
|
|
|
def GetLags(self):
|
|
|
'''
|
|
|
Function GetMoments()
|
|
|
|
|
|
Input:
|
|
|
self.dataOut.data_pre
|
|
|
self.dataOut.abscissaList
|
|
|
self.dataOut.noise
|
|
|
self.dataOut.normFactor
|
|
|
self.dataOut.data_SNR
|
|
|
self.dataOut.groupList
|
|
|
self.dataOut.nChannels
|
|
|
|
|
|
Affected:
|
|
|
self.dataOut.data_param
|
|
|
|
|
|
'''
|
|
|
|
|
|
data = self.dataOut.data_pre
|
|
|
normFactor = self.dataOut.normFactor
|
|
|
nHeights = self.dataOut.nHeights
|
|
|
absc = self.dataOut.abscissaList[:-1]
|
|
|
noise = self.dataOut.noise
|
|
|
SNR = self.dataOut.data_SNR
|
|
|
pairsList = self.dataOut.groupList
|
|
|
nChannels = self.dataOut.nChannels
|
|
|
pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
|
|
|
self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
|
|
|
|
|
|
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
|
|
|
|
|
|
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()
|
|
|
self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
|
|
|
|
|
|
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]
|
|
|
|
|
|
if (h>0)and(error1[3]<5):
|
|
|
p0 = self.dataOut.data_param[i,:,h-1]
|
|
|
else:
|
|
|
p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
|
|
|
|
|
|
try:
|
|
|
#Least Squares
|
|
|
minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
|
|
|
# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
|
|
|
#Chi square error
|
|
|
error0 = numpy.sum(infodict['fvec']**2)/(2*N)
|
|
|
#Error with Jacobian
|
|
|
error1 = self.dataOut.library.errorFunction(minp,constants,LT)
|
|
|
except:
|
|
|
minp = p0*numpy.nan
|
|
|
error0 = numpy.nan
|
|
|
error1 = p0*numpy.nan
|
|
|
|
|
|
#Save
|
|
|
if self.dataOut.data_param == None:
|
|
|
self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
|
|
|
self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
|
|
|
|
|
|
self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
|
|
|
self.dataOut.data_param[i,:,h] = minp
|
|
|
return
|
|
|
|
|
|
|
|
|
def __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
|
|
|
|
|
|
|
|
|
|
|
|
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, outputInterval):
|
|
|
|
|
|
dataTime = currentTime + paramInterval
|
|
|
deltaTime = dataTime - self.__initime
|
|
|
|
|
|
if deltaTime >= outputInterval or deltaTime < 0:
|
|
|
self.__dataReady = True
|
|
|
return
|
|
|
|
|
|
def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
|
|
|
'''
|
|
|
Function that implements winds estimation technique with detected meteors.
|
|
|
|
|
|
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.abscissaList != None:
|
|
|
absc = dataOut.abscissaList[:-1]
|
|
|
noise = dataOut.noise
|
|
|
heightList = dataOut.getHeiRange()
|
|
|
SNR = dataOut.data_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.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
|
|
|
|
|
|
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.groupList
|
|
|
nChannels = dataOut.nChannels
|
|
|
|
|
|
dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
|
|
|
dataOut.utctimeInit = dataOut.utctime
|
|
|
dataOut.outputInterval = dataOut.timeInterval
|
|
|
|
|
|
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.outputInterval = nHours*3600
|
|
|
|
|
|
if self.__isConfig == False:
|
|
|
# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
|
|
|
#Get Initial LTC time
|
|
|
self.__initime = datetime.datetime.utcfromtimestamp(self.dataOut.utctime)
|
|
|
self.__initime = (self.__initime.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.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
|
|
|
|
|
|
if self.__dataReady:
|
|
|
dataOut.utctimeInit = self.__initime
|
|
|
|
|
|
self.__initime += dataOut.outputInterval #to erase time offset
|
|
|
|
|
|
dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
|
|
|
dataOut.flagNoData = False
|
|
|
self.__buffer = None
|
|
|
|
|
|
return
|
|
|
|
|
|
class EWDriftsEstimation(Operation):
|
|
|
|
|
|
|
|
|
def __init__(self):
|
|
|
Operation.__init__(self)
|
|
|
|
|
|
def __correctValues(self, heiRang, phi, velRadial, SNR):
|
|
|
listPhi = phi.tolist()
|
|
|
maxid = listPhi.index(max(listPhi))
|
|
|
minid = listPhi.index(min(listPhi))
|
|
|
|
|
|
rango = range(len(phi))
|
|
|
# rango = numpy.delete(rango,maxid)
|
|
|
|
|
|
heiRang1 = heiRang*math.cos(phi[maxid])
|
|
|
heiRangAux = heiRang*math.cos(phi[minid])
|
|
|
indOut = (heiRang1 < heiRangAux[0]).nonzero()
|
|
|
heiRang1 = numpy.delete(heiRang1,indOut)
|
|
|
|
|
|
velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
|
|
|
SNR1 = numpy.zeros([len(phi),len(heiRang1)])
|
|
|
|
|
|
for i in rango:
|
|
|
x = heiRang*math.cos(phi[i])
|
|
|
y1 = velRadial[i,:]
|
|
|
f1 = interpolate.interp1d(x,y1,kind = 'cubic')
|
|
|
|
|
|
x1 = heiRang1
|
|
|
y11 = f1(x1)
|
|
|
|
|
|
y2 = SNR[i,:]
|
|
|
f2 = interpolate.interp1d(x,y2,kind = 'cubic')
|
|
|
y21 = f2(x1)
|
|
|
|
|
|
velRadial1[i,:] = y11
|
|
|
SNR1[i,:] = y21
|
|
|
|
|
|
return heiRang1, velRadial1, SNR1
|
|
|
|
|
|
def run(self, dataOut, zenith, zenithCorrection):
|
|
|
heiRang = dataOut.heightList
|
|
|
velRadial = dataOut.data_param[:,3,:]
|
|
|
SNR = dataOut.data_SNR
|
|
|
|
|
|
zenith = numpy.array(zenith)
|
|
|
zenith -= zenithCorrection
|
|
|
zenith *= numpy.pi/180
|
|
|
|
|
|
heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
|
|
|
|
|
|
alp = zenith[0]
|
|
|
bet = zenith[1]
|
|
|
|
|
|
w_w = velRadial1[0,:]
|
|
|
w_e = velRadial1[1,:]
|
|
|
|
|
|
w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
|
|
|
u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
|
|
|
|
|
|
winds = numpy.vstack((u,w))
|
|
|
|
|
|
dataOut.heightList = heiRang1
|
|
|
dataOut.data_output = winds
|
|
|
dataOut.data_SNR = SNR1
|
|
|
|
|
|
dataOut.utctimeInit = dataOut.utctime
|
|
|
dataOut.outputInterval = dataOut.timeInterval
|
|
|
return
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
