jrodata.py
1229 lines
| 28.4 KiB
| text/x-python
|
PythonLexer
|
r487 | ''' | ||
$Author: murco $ | ||||
$Id: JROData.py 173 2012-11-20 15:06:21Z murco $ | ||||
''' | ||||
import copy | ||||
import numpy | ||||
import datetime | ||||
from jroheaderIO import SystemHeader, RadarControllerHeader | ||||
|
r878 | from schainpy import cSchain | ||
|
r487 | |||
def getNumpyDtype(dataTypeCode): | ||||
r889 | ||||
|
r487 | if dataTypeCode == 0: | ||
numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')]) | ||||
elif dataTypeCode == 1: | ||||
numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')]) | ||||
elif dataTypeCode == 2: | ||||
numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')]) | ||||
elif dataTypeCode == 3: | ||||
numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')]) | ||||
elif dataTypeCode == 4: | ||||
numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')]) | ||||
elif dataTypeCode == 5: | ||||
numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')]) | ||||
else: | ||||
raise ValueError, 'dataTypeCode was not defined' | ||||
r889 | ||||
|
r487 | return numpyDtype | ||
def getDataTypeCode(numpyDtype): | ||||
r889 | ||||
|
r487 | if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]): | ||
datatype = 0 | ||||
elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]): | ||||
datatype = 1 | ||||
elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]): | ||||
datatype = 2 | ||||
elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]): | ||||
datatype = 3 | ||||
elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]): | ||||
datatype = 4 | ||||
elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]): | ||||
datatype = 5 | ||||
else: | ||||
datatype = None | ||||
r889 | ||||
|
r487 | return datatype | ||
def hildebrand_sekhon(data, navg): | ||||
|
r568 | """ | ||
r889 | This method is for the objective determination of the noise level in Doppler spectra. This | |||
implementation technique is based on the fact that the standard deviation of the spectral | ||||
|
r568 | densities is equal to the mean spectral density for white Gaussian noise | ||
r889 | ||||
|
r568 | Inputs: | ||
Data : heights | ||||
navg : numbers of averages | ||||
r889 | ||||
|
r568 | Return: | ||
-1 : any error | ||||
anoise : noise's level | ||||
""" | ||||
r889 | ||||
|
r931 | sortdata = numpy.sort(data, axis=None) | ||
|
r878 | # lenOfData = len(sortdata) | ||
# nums_min = lenOfData*0.2 | ||||
r889 | # | |||
|
r878 | # if nums_min <= 5: | ||
# nums_min = 5 | ||||
r889 | # | |||
|
r878 | # sump = 0. | ||
r889 | # | |||
|
r878 | # sumq = 0. | ||
r889 | # | |||
|
r878 | # j = 0 | ||
r889 | # | |||
|
r878 | # cont = 1 | ||
r889 | # | |||
|
r878 | # while((cont==1)and(j<lenOfData)): | ||
r889 | # | |||
|
r878 | # sump += sortdata[j] | ||
r889 | # | |||
|
r878 | # sumq += sortdata[j]**2 | ||
r889 | # | |||
|
r878 | # if j > nums_min: | ||
# rtest = float(j)/(j-1) + 1.0/navg | ||||
# if ((sumq*j) > (rtest*sump**2)): | ||||
# j = j - 1 | ||||
# sump = sump - sortdata[j] | ||||
# sumq = sumq - sortdata[j]**2 | ||||
# cont = 0 | ||||
r889 | # | |||
|
r878 | # j += 1 | ||
r889 | # | |||
|
r878 | # lnoise = sump /j | ||
r889 | # | |||
|
r878 | # return lnoise | ||
return cSchain.hildebrand_sekhon(sortdata, navg) | ||||
|
r487 | |||
|
r501 | class Beam: | ||
r889 | ||||
|
r501 | def __init__(self): | ||
self.codeList = [] | ||||
self.azimuthList = [] | ||||
r889 | self.zenithList = [] | |||
|
r501 | |||
|
r487 | class GenericData(object): | ||
r889 | ||||
|
r487 | flagNoData = True | ||
r889 | ||||
|
r487 | def copy(self, inputObj=None): | ||
r889 | ||||
|
r487 | if inputObj == None: | ||
return copy.deepcopy(self) | ||||
for key in inputObj.__dict__.keys(): | ||||
r889 | ||||
|
r757 | attribute = inputObj.__dict__[key] | ||
r889 | ||||
|
r757 | #If this attribute is a tuple or list | ||
if type(inputObj.__dict__[key]) in (tuple, list): | ||||
self.__dict__[key] = attribute[:] | ||||
continue | ||||
r889 | ||||
|
r757 | #If this attribute is another object or instance | ||
if hasattr(attribute, '__dict__'): | ||||
self.__dict__[key] = attribute.copy() | ||||
continue | ||||
r889 | ||||
|
r487 | self.__dict__[key] = inputObj.__dict__[key] | ||
def deepcopy(self): | ||||
r889 | ||||
|
r487 | return copy.deepcopy(self) | ||
r889 | ||||
|
r487 | def isEmpty(self): | ||
r889 | ||||
|
r487 | return self.flagNoData | ||
r889 | ||||
|
r487 | class JROData(GenericData): | ||
r889 | ||||
|
r487 | # m_BasicHeader = BasicHeader() | ||
# m_ProcessingHeader = ProcessingHeader() | ||||
systemHeaderObj = SystemHeader() | ||||
r889 | ||||
|
r487 | radarControllerHeaderObj = RadarControllerHeader() | ||
# data = None | ||||
r889 | ||||
|
r487 | type = None | ||
r889 | ||||
|
r487 | datatype = None #dtype but in string | ||
r889 | ||||
|
r487 | # dtype = None | ||
r889 | ||||
|
r487 | # nChannels = None | ||
r889 | ||||
|
r487 | # nHeights = None | ||
r889 | ||||
|
r487 | nProfiles = None | ||
r889 | ||||
|
r487 | heightList = None | ||
r889 | ||||
|
r487 | channelList = None | ||
r889 | ||||
|
r568 | flagDiscontinuousBlock = False | ||
r889 | ||||
|
r487 | useLocalTime = False | ||
r889 | ||||
|
r487 | utctime = None | ||
r889 | ||||
|
r487 | timeZone = None | ||
r889 | ||||
|
r487 | dstFlag = None | ||
r889 | ||||
|
r487 | errorCount = None | ||
r889 | ||||
|
r487 | blocksize = None | ||
r889 | ||||
|
r568 | # nCode = None | ||
r889 | # | |||
|
r568 | # nBaud = None | ||
r889 | # | |||
|
r568 | # code = None | ||
r889 | ||||
|
r487 | flagDecodeData = False #asumo q la data no esta decodificada | ||
r889 | ||||
|
r487 | flagDeflipData = False #asumo q la data no esta sin flip | ||
r889 | ||||
|
r487 | flagShiftFFT = False | ||
r889 | ||||
|
r487 | # ippSeconds = None | ||
r889 | ||||
|
r526 | # timeInterval = None | ||
r889 | ||||
|
r487 | nCohInt = None | ||
r889 | ||||
|
r568 | # noise = None | ||
r889 | ||||
|
r487 | windowOfFilter = 1 | ||
r889 | ||||
|
r487 | #Speed of ligth | ||
C = 3e8 | ||||
r889 | ||||
|
r487 | frequency = 49.92e6 | ||
r889 | ||||
|
r487 | realtime = False | ||
r889 | ||||
|
r487 | beacon_heiIndexList = None | ||
r889 | ||||
|
r487 | last_block = None | ||
r889 | ||||
|
r487 | blocknow = None | ||
|
r499 | azimuth = None | ||
r889 | ||||
|
r499 | zenith = None | ||
r889 | ||||
|
r501 | beam = Beam() | ||
r889 | ||||
|
r534 | profileIndex = None | ||
r889 | ||||
|
r487 | def getNoise(self): | ||
r889 | ||||
|
r684 | raise NotImplementedError | ||
r889 | ||||
|
r487 | def getNChannels(self): | ||
r889 | ||||
|
r487 | return len(self.channelList) | ||
r889 | ||||
|
r487 | def getChannelIndexList(self): | ||
r889 | ||||
|
r487 | return range(self.nChannels) | ||
r889 | ||||
|
r487 | def getNHeights(self): | ||
r889 | ||||
|
r487 | return len(self.heightList) | ||
r889 | ||||
|
r487 | def getHeiRange(self, extrapoints=0): | ||
r889 | ||||
|
r487 | heis = self.heightList | ||
# deltah = self.heightList[1] - self.heightList[0] | ||||
r889 | # | |||
|
r487 | # heis.append(self.heightList[-1]) | ||
r889 | ||||
|
r487 | return heis | ||
r889 | ||||
|
r765 | def getDeltaH(self): | ||
r889 | ||||
|
r765 | delta = self.heightList[1] - self.heightList[0] | ||
r889 | ||||
|
r765 | return delta | ||
r889 | ||||
|
r487 | def getltctime(self): | ||
r889 | ||||
|
r487 | if self.useLocalTime: | ||
return self.utctime - self.timeZone*60 | ||||
r889 | ||||
|
r487 | return self.utctime | ||
r889 | ||||
|
r487 | def getDatatime(self): | ||
r889 | ||||
|
r487 | datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime) | ||
return datatimeValue | ||||
r889 | ||||
|
r487 | def getTimeRange(self): | ||
r889 | ||||
|
r487 | datatime = [] | ||
r889 | ||||
|
r487 | datatime.append(self.ltctime) | ||
|
r815 | datatime.append(self.ltctime + self.timeInterval+1) | ||
r889 | ||||
|
r487 | datatime = numpy.array(datatime) | ||
r889 | ||||
|
r487 | return datatime | ||
r889 | ||||
|
r765 | def getFmaxTimeResponse(self): | ||
r889 | ||||
|
r765 | period = (10**-6)*self.getDeltaH()/(0.15) | ||
r889 | ||||
|
r765 | PRF = 1./(period * self.nCohInt) | ||
r889 | ||||
|
r765 | fmax = PRF | ||
r889 | ||||
|
r765 | return fmax | ||
r889 | ||||
|
r487 | def getFmax(self): | ||
r889 | ||||
|
r487 | PRF = 1./(self.ippSeconds * self.nCohInt) | ||
r889 | ||||
|
r765 | fmax = PRF | ||
r889 | ||||
|
r487 | return fmax | ||
r889 | ||||
|
r487 | def getVmax(self): | ||
r889 | ||||
|
r487 | _lambda = self.C/self.frequency | ||
r889 | ||||
|
r850 | vmax = self.getFmax() * _lambda/2 | ||
r889 | ||||
|
r487 | return vmax | ||
r889 | ||||
|
r487 | def get_ippSeconds(self): | ||
''' | ||||
''' | ||||
return self.radarControllerHeaderObj.ippSeconds | ||||
r889 | ||||
|
r487 | def set_ippSeconds(self, ippSeconds): | ||
''' | ||||
''' | ||||
r889 | ||||
|
r487 | self.radarControllerHeaderObj.ippSeconds = ippSeconds | ||
r889 | ||||
|
r487 | return | ||
r889 | ||||
|
r487 | def get_dtype(self): | ||
''' | ||||
''' | ||||
return getNumpyDtype(self.datatype) | ||||
r889 | ||||
|
r487 | def set_dtype(self, numpyDtype): | ||
''' | ||||
''' | ||||
r889 | ||||
|
r487 | self.datatype = getDataTypeCode(numpyDtype) | ||
|
r568 | |||
def get_code(self): | ||||
''' | ||||
''' | ||||
return self.radarControllerHeaderObj.code | ||||
r889 | ||||
|
r568 | def set_code(self, code): | ||
''' | ||||
''' | ||||
self.radarControllerHeaderObj.code = code | ||||
r889 | ||||
|
r568 | return | ||
def get_ncode(self): | ||||
''' | ||||
''' | ||||
return self.radarControllerHeaderObj.nCode | ||||
r889 | ||||
|
r568 | def set_ncode(self, nCode): | ||
''' | ||||
''' | ||||
self.radarControllerHeaderObj.nCode = nCode | ||||
r889 | ||||
|
r568 | return | ||
def get_nbaud(self): | ||||
''' | ||||
''' | ||||
return self.radarControllerHeaderObj.nBaud | ||||
r889 | ||||
|
r568 | def set_nbaud(self, nBaud): | ||
''' | ||||
''' | ||||
self.radarControllerHeaderObj.nBaud = nBaud | ||||
r889 | ||||
|
r568 | return | ||
r889 | ||||
|
r487 | nChannels = property(getNChannels, "I'm the 'nChannel' property.") | ||
channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.") | ||||
nHeights = property(getNHeights, "I'm the 'nHeights' property.") | ||||
#noise = property(getNoise, "I'm the 'nHeights' property.") | ||||
datatime = property(getDatatime, "I'm the 'datatime' property") | ||||
ltctime = property(getltctime, "I'm the 'ltctime' property") | ||||
ippSeconds = property(get_ippSeconds, set_ippSeconds) | ||||
dtype = property(get_dtype, set_dtype) | ||||
|
r526 | # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property") | ||
|
r568 | code = property(get_code, set_code) | ||
nCode = property(get_ncode, set_ncode) | ||||
nBaud = property(get_nbaud, set_nbaud) | ||||
r889 | ||||
|
r487 | class Voltage(JROData): | ||
r889 | ||||
|
r487 | #data es un numpy array de 2 dmensiones (canales, alturas) | ||
data = None | ||||
r889 | ||||
|
r487 | def __init__(self): | ||
''' | ||||
Constructor | ||||
''' | ||||
r889 | ||||
|
r566 | self.useLocalTime = True | ||
r889 | ||||
|
r487 | self.radarControllerHeaderObj = RadarControllerHeader() | ||
r889 | ||||
|
r487 | self.systemHeaderObj = SystemHeader() | ||
r889 | ||||
|
r487 | self.type = "Voltage" | ||
r889 | ||||
|
r487 | self.data = None | ||
r889 | ||||
|
r487 | # self.dtype = None | ||
r889 | ||||
|
r487 | # self.nChannels = 0 | ||
r889 | ||||
|
r487 | # self.nHeights = 0 | ||
r889 | ||||
|
r487 | self.nProfiles = None | ||
r889 | ||||
|
r487 | self.heightList = None | ||
r889 | ||||
|
r487 | self.channelList = None | ||
r889 | ||||
|
r487 | # self.channelIndexList = None | ||
r889 | ||||
|
r487 | self.flagNoData = True | ||
r889 | ||||
|
r568 | self.flagDiscontinuousBlock = False | ||
r889 | ||||
|
r487 | self.utctime = None | ||
r889 | ||||
|
r487 | self.timeZone = None | ||
r889 | ||||
|
r487 | self.dstFlag = None | ||
r889 | ||||
|
r487 | self.errorCount = None | ||
r889 | ||||
|
r487 | self.nCohInt = None | ||
r889 | ||||
|
r487 | self.blocksize = None | ||
r889 | ||||
|
r487 | self.flagDecodeData = False #asumo q la data no esta decodificada | ||
r889 | ||||
|
r487 | self.flagDeflipData = False #asumo q la data no esta sin flip | ||
r889 | ||||
|
r487 | self.flagShiftFFT = False | ||
r889 | ||||
|
r526 | self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil | ||
r889 | ||||
|
r534 | self.profileIndex = 0 | ||
r889 | ||||
|
r568 | def getNoisebyHildebrand(self, channel = None): | ||
|
r487 | """ | ||
Determino el nivel de ruido usando el metodo Hildebrand-Sekhon | ||||
r889 | ||||
|
r487 | Return: | ||
noiselevel | ||||
""" | ||||
|
r568 | if channel != None: | ||
data = self.data[channel] | ||||
nChannels = 1 | ||||
else: | ||||
data = self.data | ||||
nChannels = self.nChannels | ||||
r889 | ||||
|
r568 | noise = numpy.zeros(nChannels) | ||
power = data * numpy.conjugate(data) | ||||
r889 | ||||
|
r568 | for thisChannel in range(nChannels): | ||
if nChannels == 1: | ||||
daux = power[:].real | ||||
else: | ||||
daux = power[thisChannel,:].real | ||||
noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt) | ||||
r889 | ||||
|
r568 | return noise | ||
r889 | ||||
|
r568 | def getNoise(self, type = 1, channel = None): | ||
r889 | ||||
|
r487 | if type == 1: | ||
|
r568 | noise = self.getNoisebyHildebrand(channel) | ||
r889 | ||||
|
r720 | return noise | ||
r889 | ||||
|
r568 | def getPower(self, channel = None): | ||
r889 | ||||
|
r568 | if channel != None: | ||
data = self.data[channel] | ||||
else: | ||||
data = self.data | ||||
r889 | ||||
|
r568 | power = data * numpy.conjugate(data) | ||
|
r850 | powerdB = 10*numpy.log10(power.real) | ||
powerdB = numpy.squeeze(powerdB) | ||||
r889 | ||||
|
r850 | return powerdB | ||
r889 | ||||
|
r526 | def getTimeInterval(self): | ||
r889 | ||||
|
r526 | timeInterval = self.ippSeconds * self.nCohInt | ||
r889 | ||||
|
r526 | return timeInterval | ||
r889 | ||||
|
r526 | noise = property(getNoise, "I'm the 'nHeights' property.") | ||
timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property") | ||||
r889 | ||||
|
r487 | class Spectra(JROData): | ||
r889 | ||||
|
r832 | #data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas) | ||
|
r487 | data_spc = None | ||
r889 | ||||
|
r832 | #data cspc es un numpy array de 2 dmensiones (canales, pares, alturas) | ||
|
r487 | data_cspc = None | ||
r889 | ||||
|
r832 | #data dc es un numpy array de 2 dmensiones (canales, alturas) | ||
|
r487 | data_dc = None | ||
r889 | ||||
|
r832 | #data power | ||
data_pwr = None | ||||
r889 | ||||
|
r487 | nFFTPoints = None | ||
r889 | ||||
|
r487 | # nPairs = None | ||
r889 | ||||
|
r487 | pairsList = None | ||
r889 | ||||
|
r487 | nIncohInt = None | ||
r889 | ||||
|
r487 | wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia | ||
r889 | ||||
|
r487 | nCohInt = None #se requiere para determinar el valor de timeInterval | ||
r889 | ||||
|
r487 | ippFactor = None | ||
r889 | ||||
|
r534 | profileIndex = 0 | ||
r889 | ||||
|
r773 | plotting = "spectra" | ||
r889 | ||||
|
r487 | def __init__(self): | ||
''' | ||||
Constructor | ||||
''' | ||||
r889 | ||||
|
r566 | self.useLocalTime = True | ||
r889 | ||||
|
r487 | self.radarControllerHeaderObj = RadarControllerHeader() | ||
r889 | ||||
|
r487 | self.systemHeaderObj = SystemHeader() | ||
r889 | ||||
|
r487 | self.type = "Spectra" | ||
r889 | ||||
|
r487 | # self.data = None | ||
r889 | ||||
|
r487 | # self.dtype = None | ||
r889 | ||||
|
r487 | # self.nChannels = 0 | ||
r889 | ||||
|
r487 | # self.nHeights = 0 | ||
r889 | ||||
|
r487 | self.nProfiles = None | ||
r889 | ||||
|
r487 | self.heightList = None | ||
r889 | ||||
|
r487 | self.channelList = None | ||
r889 | ||||
|
r487 | # self.channelIndexList = None | ||
self.pairsList = None | ||||
r889 | ||||
|
r487 | self.flagNoData = True | ||
r889 | ||||
|
r568 | self.flagDiscontinuousBlock = False | ||
r889 | ||||
|
r487 | self.utctime = None | ||
r889 | ||||
|
r487 | self.nCohInt = None | ||
r889 | ||||
|
r487 | self.nIncohInt = None | ||
r889 | ||||
|
r487 | self.blocksize = None | ||
r889 | ||||
|
r487 | self.nFFTPoints = None | ||
r889 | ||||
|
r487 | self.wavelength = None | ||
r889 | ||||
|
r487 | self.flagDecodeData = False #asumo q la data no esta decodificada | ||
r889 | ||||
|
r487 | self.flagDeflipData = False #asumo q la data no esta sin flip | ||
r889 | ||||
|
r487 | self.flagShiftFFT = False | ||
r889 | ||||
|
r487 | self.ippFactor = 1 | ||
r889 | ||||
|
r487 | #self.noise = None | ||
r889 | ||||
|
r487 | self.beacon_heiIndexList = [] | ||
r889 | ||||
|
r487 | self.noise_estimation = None | ||
r889 | ||||
|
r568 | def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None): | ||
|
r487 | """ | ||
Determino el nivel de ruido usando el metodo Hildebrand-Sekhon | ||||
r889 | ||||
|
r487 | Return: | ||
noiselevel | ||||
""" | ||||
r889 | ||||
|
r487 | noise = numpy.zeros(self.nChannels) | ||
r889 | ||||
|
r487 | for channel in range(self.nChannels): | ||
|
r568 | daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index] | ||
|
r487 | noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) | ||
r889 | ||||
return noise | ||||
|
r568 | def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None): | ||
r889 | ||||
|
r729 | if self.noise_estimation is not None: | ||
|
r487 | return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py | ||
else: | ||||
|
r568 | noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index) | ||
|
r487 | return noise | ||
r889 | ||||
|
r765 | def getFreqRangeTimeResponse(self, extrapoints=0): | ||
r889 | ||||
|
r765 | deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints*self.ippFactor) | ||
freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2 | ||||
r889 | ||||
|
r771 | return freqrange | ||
r889 | ||||
|
r771 | def getAcfRange(self, extrapoints=0): | ||
r889 | ||||
|
r771 | deltafreq = 10./(self.getFmax() / (self.nFFTPoints*self.ippFactor)) | ||
freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2 | ||||
r889 | ||||
|
r765 | return freqrange | ||
r889 | ||||
|
r487 | def getFreqRange(self, extrapoints=0): | ||
r889 | ||||
|
r487 | deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor) | ||
freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2 | ||||
r889 | ||||
|
r487 | return freqrange | ||
def getVelRange(self, extrapoints=0): | ||||
|
r1146 | |||
|
r487 | deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor) | ||
r889 | velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) #- deltav/2 | |||
|
r487 | return velrange | ||
r889 | ||||
|
r487 | def getNPairs(self): | ||
r889 | ||||
|
r487 | return len(self.pairsList) | ||
r889 | ||||
|
r487 | def getPairsIndexList(self): | ||
r889 | ||||
|
r487 | return range(self.nPairs) | ||
r889 | ||||
|
r487 | def getNormFactor(self): | ||
r889 | ||||
|
r487 | pwcode = 1 | ||
r889 | ||||
|
r487 | if self.flagDecodeData: | ||
pwcode = numpy.sum(self.code[0]**2) | ||||
#normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter | ||||
normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter | ||||
r889 | ||||
|
r487 | return normFactor | ||
r889 | ||||
|
r487 | def getFlagCspc(self): | ||
r889 | ||||
|
r611 | if self.data_cspc is None: | ||
|
r487 | return True | ||
r889 | ||||
|
r487 | return False | ||
r889 | ||||
|
r487 | def getFlagDc(self): | ||
r889 | ||||
|
r611 | if self.data_dc is None: | ||
|
r487 | return True | ||
r889 | ||||
|
r487 | return False | ||
r889 | ||||
|
r526 | def getTimeInterval(self): | ||
r889 | ||||
|
r526 | timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles | ||
r889 | ||||
|
r526 | return timeInterval | ||
r889 | ||||
|
r832 | def getPower(self): | ||
r889 | ||||
|
r832 | factor = self.normFactor | ||
z = self.data_spc/factor | ||||
r889 | z = numpy.where(numpy.isfinite(z), z, numpy.NAN) | |||
|
r832 | avg = numpy.average(z, axis=1) | ||
r889 | ||||
|
r832 | return 10*numpy.log10(avg) | ||
r889 | ||||
def getCoherence(self, pairsList=None, phase=False): | ||||
z = [] | ||||
if pairsList is None: | ||||
pairsIndexList = self.pairsIndexList | ||||
else: | ||||
pairsIndexList = [] | ||||
for pair in pairsList: | ||||
if pair not in self.pairsList: | ||||
raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair) | ||||
pairsIndexList.append(self.pairsList.index(pair)) | ||||
for i in range(len(pairsIndexList)): | ||||
pair = self.pairsList[pairsIndexList[i]] | ||||
ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0) | ||||
powa = numpy.average(self.data_spc[pair[0], :, :], axis=0) | ||||
powb = numpy.average(self.data_spc[pair[1], :, :], axis=0) | ||||
avgcoherenceComplex = ccf/numpy.sqrt(powa*powb) | ||||
if phase: | ||||
data = numpy.arctan2(avgcoherenceComplex.imag, | ||||
avgcoherenceComplex.real)*180/numpy.pi | ||||
else: | ||||
data = numpy.abs(avgcoherenceComplex) | ||||
z.append(data) | ||||
return numpy.array(z) | ||||
|
r624 | def setValue(self, value): | ||
r889 | ||||
|
r624 | print "This property should not be initialized" | ||
r889 | ||||
|
r624 | return | ||
r889 | ||||
|
r624 | nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.") | ||
pairsIndexList = property(getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.") | ||||
normFactor = property(getNormFactor, setValue, "I'm the 'getNormFactor' property.") | ||||
flag_cspc = property(getFlagCspc, setValue) | ||||
flag_dc = property(getFlagDc, setValue) | ||||
noise = property(getNoise, setValue, "I'm the 'nHeights' property.") | ||||
timeInterval = property(getTimeInterval, setValue, "I'm the 'timeInterval' property") | ||||
r889 | ||||
|
r487 | class SpectraHeis(Spectra): | ||
r889 | ||||
|
r487 | data_spc = None | ||
r889 | ||||
|
r487 | data_cspc = None | ||
r889 | ||||
|
r487 | data_dc = None | ||
r889 | ||||
|
r487 | nFFTPoints = None | ||
r889 | ||||
|
r487 | # nPairs = None | ||
r889 | ||||
|
r487 | pairsList = None | ||
r889 | ||||
|
r587 | nCohInt = None | ||
r889 | ||||
|
r487 | nIncohInt = None | ||
r889 | ||||
|
r487 | def __init__(self): | ||
r889 | ||||
|
r487 | self.radarControllerHeaderObj = RadarControllerHeader() | ||
r889 | ||||
|
r487 | self.systemHeaderObj = SystemHeader() | ||
r889 | ||||
|
r487 | self.type = "SpectraHeis" | ||
r889 | ||||
|
r487 | # self.dtype = None | ||
r889 | ||||
|
r487 | # self.nChannels = 0 | ||
r889 | ||||
|
r487 | # self.nHeights = 0 | ||
r889 | ||||
|
r487 | self.nProfiles = None | ||
r889 | ||||
|
r487 | self.heightList = None | ||
r889 | ||||
|
r487 | self.channelList = None | ||
r889 | ||||
|
r487 | # self.channelIndexList = None | ||
r889 | ||||
|
r487 | self.flagNoData = True | ||
r889 | ||||
|
r568 | self.flagDiscontinuousBlock = False | ||
r889 | ||||
|
r487 | # self.nPairs = 0 | ||
r889 | ||||
|
r487 | self.utctime = None | ||
r889 | ||||
|
r487 | self.blocksize = None | ||
r889 | ||||
|
r534 | self.profileIndex = 0 | ||
r889 | ||||
|
r587 | self.nCohInt = 1 | ||
r889 | ||||
|
r587 | self.nIncohInt = 1 | ||
r889 | ||||
|
r496 | def getNormFactor(self): | ||
pwcode = 1 | ||||
if self.flagDecodeData: | ||||
pwcode = numpy.sum(self.code[0]**2) | ||||
r889 | ||||
|
r496 | normFactor = self.nIncohInt*self.nCohInt*pwcode | ||
r889 | ||||
|
r496 | return normFactor | ||
r889 | ||||
|
r526 | def getTimeInterval(self): | ||
r889 | ||||
|
r526 | timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt | ||
r889 | ||||
|
r526 | return timeInterval | ||
r889 | ||||
|
r496 | normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.") | ||
|
r526 | timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property") | ||
|
r487 | |||
|
r587 | class Fits(JROData): | ||
r889 | ||||
|
r487 | heightList = None | ||
r889 | ||||
|
r487 | channelList = None | ||
r889 | ||||
|
r487 | flagNoData = True | ||
r889 | ||||
|
r568 | flagDiscontinuousBlock = False | ||
r889 | ||||
|
r487 | useLocalTime = False | ||
r889 | ||||
|
r487 | utctime = None | ||
r889 | ||||
|
r487 | timeZone = None | ||
r889 | ||||
|
r487 | # ippSeconds = None | ||
r889 | ||||
|
r526 | # timeInterval = None | ||
r889 | ||||
|
r487 | nCohInt = None | ||
r889 | ||||
|
r487 | nIncohInt = None | ||
r889 | ||||
|
r487 | noise = None | ||
r889 | ||||
|
r487 | windowOfFilter = 1 | ||
r889 | ||||
|
r487 | #Speed of ligth | ||
C = 3e8 | ||||
r889 | ||||
|
r487 | frequency = 49.92e6 | ||
r889 | ||||
|
r487 | realtime = False | ||
r889 | ||||
|
r487 | def __init__(self): | ||
r889 | ||||
|
r487 | self.type = "Fits" | ||
r889 | ||||
|
r487 | self.nProfiles = None | ||
r889 | ||||
|
r487 | self.heightList = None | ||
r889 | ||||
|
r487 | self.channelList = None | ||
r889 | ||||
|
r487 | # self.channelIndexList = None | ||
r889 | ||||
|
r487 | self.flagNoData = True | ||
r889 | ||||
|
r487 | self.utctime = None | ||
r889 | ||||
|
r587 | self.nCohInt = 1 | ||
r889 | ||||
|
r587 | self.nIncohInt = 1 | ||
r889 | ||||
|
r487 | self.useLocalTime = True | ||
r889 | ||||
|
r534 | self.profileIndex = 0 | ||
r889 | ||||
|
r487 | # self.utctime = None | ||
# self.timeZone = None | ||||
# self.ltctime = None | ||||
# self.timeInterval = None | ||||
# self.header = None | ||||
# self.data_header = None | ||||
# self.data = None | ||||
# self.datatime = None | ||||
# self.flagNoData = False | ||||
# self.expName = '' | ||||
# self.nChannels = None | ||||
# self.nSamples = None | ||||
# self.dataBlocksPerFile = None | ||||
# self.comments = '' | ||||
r889 | # | |||
|
r487 | |||
def getltctime(self): | ||||
r889 | ||||
|
r487 | if self.useLocalTime: | ||
return self.utctime - self.timeZone*60 | ||||
r889 | ||||
|
r487 | return self.utctime | ||
r889 | ||||
|
r487 | def getDatatime(self): | ||
r889 | ||||
|
r487 | datatime = datetime.datetime.utcfromtimestamp(self.ltctime) | ||
return datatime | ||||
r889 | ||||
|
r487 | def getTimeRange(self): | ||
r889 | ||||
|
r487 | datatime = [] | ||
r889 | ||||
|
r487 | datatime.append(self.ltctime) | ||
datatime.append(self.ltctime + self.timeInterval) | ||||
r889 | ||||
|
r487 | datatime = numpy.array(datatime) | ||
r889 | ||||
|
r487 | return datatime | ||
r889 | ||||
|
r487 | def getHeiRange(self): | ||
r889 | ||||
|
r487 | heis = self.heightList | ||
r889 | ||||
|
r487 | return heis | ||
r889 | ||||
|
r487 | def getNHeights(self): | ||
r889 | ||||
|
r487 | return len(self.heightList) | ||
r889 | ||||
|
r487 | def getNChannels(self): | ||
r889 | ||||
|
r487 | return len(self.channelList) | ||
r889 | ||||
|
r487 | def getChannelIndexList(self): | ||
r889 | ||||
|
r487 | return range(self.nChannels) | ||
r889 | ||||
|
r487 | def getNoise(self, type = 1): | ||
r889 | ||||
|
r568 | #noise = numpy.zeros(self.nChannels) | ||
r889 | ||||
|
r487 | if type == 1: | ||
noise = self.getNoisebyHildebrand() | ||||
r889 | ||||
|
r487 | if type == 2: | ||
noise = self.getNoisebySort() | ||||
r889 | ||||
|
r487 | if type == 3: | ||
noise = self.getNoisebyWindow() | ||||
r889 | ||||
|
r487 | return noise | ||
r889 | ||||
|
r526 | def getTimeInterval(self): | ||
r889 | ||||
|
r587 | timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt | ||
r889 | ||||
|
r587 | return timeInterval | ||
r889 | ||||
|
r487 | datatime = property(getDatatime, "I'm the 'datatime' property") | ||
nHeights = property(getNHeights, "I'm the 'nHeights' property.") | ||||
nChannels = property(getNChannels, "I'm the 'nChannel' property.") | ||||
channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.") | ||||
noise = property(getNoise, "I'm the 'nHeights' property.") | ||||
r889 | ||||
|
r487 | ltctime = property(getltctime, "I'm the 'ltctime' property") | ||
|
r526 | timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property") | ||
r889 | ||||
|
r502 | class Correlation(JROData): | ||
r889 | ||||
|
r502 | noise = None | ||
r889 | ||||
|
r502 | SNR = None | ||
r889 | ||||
|
r502 | #-------------------------------------------------- | ||
r889 | ||||
|
r850 | mode = None | ||
r889 | ||||
|
r850 | split = False | ||
data_cf = None | ||||
r889 | ||||
|
r850 | lags = None | ||
r889 | ||||
|
r850 | lagRange = None | ||
r889 | ||||
|
r502 | pairsList = None | ||
r889 | ||||
|
r850 | normFactor = None | ||
r889 | ||||
|
r850 | #-------------------------------------------------- | ||
r889 | ||||
|
r850 | # calculateVelocity = None | ||
r889 | ||||
|
r850 | nLags = None | ||
r889 | ||||
|
r850 | nPairs = None | ||
r889 | ||||
|
r850 | nAvg = None | ||
r889 | ||||
|
r502 | def __init__(self): | ||
''' | ||||
Constructor | ||||
''' | ||||
self.radarControllerHeaderObj = RadarControllerHeader() | ||||
r889 | ||||
|
r502 | self.systemHeaderObj = SystemHeader() | ||
r889 | ||||
|
r502 | self.type = "Correlation" | ||
r889 | ||||
|
r502 | self.data = None | ||
r889 | ||||
|
r502 | self.dtype = None | ||
r889 | ||||
|
r502 | self.nProfiles = None | ||
r889 | ||||
|
r502 | self.heightList = None | ||
r889 | ||||
|
r502 | self.channelList = None | ||
r889 | ||||
|
r502 | self.flagNoData = True | ||
r889 | ||||
|
r568 | self.flagDiscontinuousBlock = False | ||
r889 | ||||
|
r502 | self.utctime = None | ||
r889 | ||||
|
r502 | self.timeZone = None | ||
r889 | ||||
|
r502 | self.dstFlag = None | ||
r889 | ||||
|
r502 | self.errorCount = None | ||
r889 | ||||
|
r502 | self.blocksize = None | ||
r889 | ||||
|
r502 | self.flagDecodeData = False #asumo q la data no esta decodificada | ||
r889 | ||||
|
r502 | self.flagDeflipData = False #asumo q la data no esta sin flip | ||
r889 | ||||
|
r502 | self.pairsList = None | ||
r889 | ||||
|
r502 | self.nPoints = None | ||
|
r850 | |||
|
r502 | def getPairsList(self): | ||
r889 | ||||
|
r502 | return self.pairsList | ||
r889 | ||||
|
r502 | def getNoise(self, mode = 2): | ||
r889 | ||||
|
r502 | indR = numpy.where(self.lagR == 0)[0][0] | ||
indT = numpy.where(self.lagT == 0)[0][0] | ||||
r889 | ||||
|
r502 | jspectra0 = self.data_corr[:,:,indR,:] | ||
jspectra = copy.copy(jspectra0) | ||||
r889 | ||||
|
r502 | num_chan = jspectra.shape[0] | ||
num_hei = jspectra.shape[2] | ||||
r889 | ||||
|
r502 | freq_dc = jspectra.shape[1]/2 | ||
r889 | ind_vel = numpy.array([-2,-1,1,2]) + freq_dc | |||
|
r502 | if ind_vel[0]<0: | ||
ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof | ||||
r889 | ||||
if mode == 1: | ||||
|
r502 | jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION | ||
r889 | ||||
|
r502 | if mode == 2: | ||
r889 | ||||
|
r502 | vel = numpy.array([-2,-1,1,2]) | ||
xx = numpy.zeros([4,4]) | ||||
r889 | ||||
|
r502 | for fil in range(4): | ||
xx[fil,:] = vel[fil]**numpy.asarray(range(4)) | ||||
r889 | ||||
|
r502 | xx_inv = numpy.linalg.inv(xx) | ||
xx_aux = xx_inv[0,:] | ||||
r889 | ||||
|
r502 | for ich in range(num_chan): | ||
yy = jspectra[ich,ind_vel,:] | ||||
jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy) | ||||
junkid = jspectra[ich,freq_dc,:]<=0 | ||||
cjunkid = sum(junkid) | ||||
r889 | ||||
|
r502 | if cjunkid.any(): | ||
jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2 | ||||
r889 | ||||
|
r502 | noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:] | ||
r889 | ||||
|
r502 | return noise | ||
|
r587 | |||
def getTimeInterval(self): | ||||
r889 | ||||
|
r850 | timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles | ||
r889 | ||||
|
r587 | return timeInterval | ||
r889 | ||||
|
r850 | def splitFunctions(self): | ||
r889 | ||||
|
r850 | pairsList = self.pairsList | ||
ccf_pairs = [] | ||||
acf_pairs = [] | ||||
ccf_ind = [] | ||||
acf_ind = [] | ||||
r889 | for l in range(len(pairsList)): | |||
|
r850 | chan0 = pairsList[l][0] | ||
chan1 = pairsList[l][1] | ||||
r889 | ||||
#Obteniendo pares de Autocorrelacion | ||||
|
r850 | if chan0 == chan1: | ||
acf_pairs.append(chan0) | ||||
acf_ind.append(l) | ||||
else: | ||||
ccf_pairs.append(pairsList[l]) | ||||
ccf_ind.append(l) | ||||
r889 | ||||
|
r850 | data_acf = self.data_cf[acf_ind] | ||
data_ccf = self.data_cf[ccf_ind] | ||||
return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf | ||||
def getNormFactor(self): | ||||
acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions() | ||||
acf_pairs = numpy.array(acf_pairs) | ||||
normFactor = numpy.zeros((self.nPairs,self.nHeights)) | ||||
r889 | ||||
|
r850 | for p in range(self.nPairs): | ||
pair = self.pairsList[p] | ||||
r889 | ||||
|
r850 | ch0 = pair[0] | ||
ch1 = pair[1] | ||||
r889 | ||||
|
r850 | ch0_max = numpy.max(data_acf[acf_pairs==ch0,:,:], axis=1) | ||
ch1_max = numpy.max(data_acf[acf_pairs==ch1,:,:], axis=1) | ||||
normFactor[p,:] = numpy.sqrt(ch0_max*ch1_max) | ||||
r889 | ||||
|
r850 | return normFactor | ||
r889 | ||||
|
r587 | timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property") | ||
|
r850 | normFactor = property(getNormFactor, "I'm the 'normFactor property'") | ||
r889 | ||||
r922 | class Parameters(Spectra): | |||
|
r502 | |||
|
r832 | experimentInfo = None #Information about the experiment | ||
r889 | ||||
|
r513 | #Information from previous data | ||
r889 | ||||
|
r502 | inputUnit = None #Type of data to be processed | ||
r889 | ||||
|
r502 | operation = None #Type of operation to parametrize | ||
r889 | ||||
r922 | #normFactor = None #Normalization Factor | |||
r889 | ||||
|
r513 | groupList = None #List of Pairs, Groups, etc | ||
r889 | ||||
|
r513 | #Parameters | ||
r889 | ||||
|
r502 | data_param = None #Parameters obtained | ||
r889 | ||||
|
r502 | data_pre = None #Data Pre Parametrization | ||
r889 | ||||
|
r514 | data_SNR = None #Signal to Noise Ratio | ||
r889 | ||||
|
r608 | # heightRange = None #Heights | ||
r889 | ||||
|
r608 | abscissaList = None #Abscissa, can be velocities, lags or time | ||
r889 | ||||
r923 | # noise = None #Noise Potency | |||
r889 | ||||
|
r608 | utctimeInit = None #Initial UTC time | ||
r889 | ||||
|
r502 | paramInterval = None #Time interval to calculate Parameters in seconds | ||
r889 | ||||
|
r802 | useLocalTime = True | ||
r889 | ||||
|
r513 | #Fitting | ||
r889 | ||||
|
r514 | data_error = None #Error of the estimation | ||
r889 | ||||
constants = None | ||||
|
r513 | library = None | ||
r889 | ||||
|
r513 | #Output signal | ||
r889 | ||||
|
r513 | outputInterval = None #Time interval to calculate output signal in seconds | ||
r889 | ||||
|
r513 | data_output = None #Out signal | ||
r889 | ||||
|
r850 | nAvg = None | ||
r889 | ||||
r922 | noise_estimation = None | |||
|
r1001 | |||
GauSPC = None #Fit gaussian SPC | ||||
r922 | ||||
r889 | ||||
|
r502 | def __init__(self): | ||
''' | ||||
Constructor | ||||
''' | ||||
self.radarControllerHeaderObj = RadarControllerHeader() | ||||
r889 | ||||
|
r502 | self.systemHeaderObj = SystemHeader() | ||
r889 | ||||
|
r502 | self.type = "Parameters" | ||
r889 | ||||
|
r804 | def getTimeRange1(self, interval): | ||
r889 | ||||
|
r502 | datatime = [] | ||
r889 | ||||
|
r608 | if self.useLocalTime: | ||
time1 = self.utctimeInit - self.timeZone*60 | ||||
else: | ||||
|
r802 | time1 = self.utctimeInit | ||
|
r1018 | |||
|
r608 | datatime.append(time1) | ||
|
r804 | datatime.append(time1 + interval) | ||
|
r502 | datatime = numpy.array(datatime) | ||
r889 | ||||
r922 | return datatime | |||
r923 | ||||
def getTimeInterval(self): | ||||
|
r937 | if hasattr(self, 'timeInterval1'): | ||
return self.timeInterval1 | ||||
else: | ||||
return self.paramInterval | ||||
|
r928 | |||
|
r1001 | def setValue(self, value): | ||
print "This property should not be initialized" | ||||
return | ||||
|
r928 | def getNoise(self): | ||
return self.spc_noise | ||||
timeInterval = property(getTimeInterval) | ||||
|
r1001 | noise = property(getNoise, setValue, "I'm the 'Noise' property.") | ||