jrodata.py
1233 lines
| 28.9 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 | |||
|
r1167 | from .jroheaderIO import SystemHeader, RadarControllerHeader | |
|
r878 | ||
|
r487 | ||
def getNumpyDtype(dataTypeCode): | |||
r889 | |||
|
r487 | if dataTypeCode == 0: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')]) | |
|
r487 | elif dataTypeCode == 1: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')]) | |
|
r487 | elif dataTypeCode == 2: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')]) | |
|
r487 | elif dataTypeCode == 3: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')]) | |
|
r487 | elif dataTypeCode == 4: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')]) | |
|
r487 | elif dataTypeCode == 5: | |
|
r1092 | numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')]) | |
|
r487 | else: | |
|
r1167 | raise ValueError('dataTypeCode was not defined') | |
r889 | |||
|
r487 | return numpyDtype | |
|
r1092 | ||
|
r487 | def getDataTypeCode(numpyDtype): | |
r889 | |||
|
r1092 | if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]): | |
|
r487 | datatype = 0 | |
|
r1092 | elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]): | |
|
r487 | datatype = 1 | |
|
r1092 | elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]): | |
|
r487 | datatype = 2 | |
|
r1092 | elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]): | |
|
r487 | datatype = 3 | |
|
r1092 | elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]): | |
|
r487 | datatype = 4 | |
|
r1092 | elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]): | |
|
r487 | datatype = 5 | |
else: | |||
datatype = None | |||
r889 | |||
|
r487 | return datatype | |
|
r1092 | ||
|
r487 | 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: | |
|
r1176 | mean : noise's level | |
|
r568 | """ | |
r889 | |||
|
r1176 | sorted_spectrum = numpy.sort(data, axis=None) | |
nnoise = len(sorted_spectrum) # default to all points in the spectrum as noise | |||
for npts in range(1, len(sorted_spectrum)+1): | |||
partial = sorted_spectrum[:npts] | |||
mean = partial.mean() | |||
var = partial.var() | |||
if var * navg < mean**2.: | |||
nnoise = npts | |||
else: | |||
# partial spectrum no longer has characteristics of white noise | |||
break | |||
|
r1167 | ||
|
r1176 | noise_spectrum = sorted_spectrum[:nnoise] | |
mean = noise_spectrum.mean() | |||
return mean | |||
|
r487 | ||
|
r501 | class Beam: | |
r889 | |||
|
r501 | def __init__(self): | |
self.codeList = [] | |||
self.azimuthList = [] | |||
r889 | self.zenithList = [] | ||
|
r501 | ||
|
r1092 | ||
|
r487 | class GenericData(object): | |
r889 | |||
|
r487 | flagNoData = True | |
r889 | |||
|
r487 | def copy(self, inputObj=None): | |
r889 | |||
|
r487 | if inputObj == None: | |
return copy.deepcopy(self) | |||
|
r1167 | for key in list(inputObj.__dict__.keys()): | |
r889 | |||
|
r757 | attribute = inputObj.__dict__[key] | |
r889 | |||
|
r1092 | # If this attribute is a tuple or list | |
|
r757 | if type(inputObj.__dict__[key]) in (tuple, list): | |
self.__dict__[key] = attribute[:] | |||
continue | |||
r889 | |||
|
r1092 | # If this attribute is another object or instance | |
|
r757 | 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 | |||
|
r1092 | ||
|
r487 | class JROData(GenericData): | |
r889 | |||
|
r1092 | # m_BasicHeader = BasicHeader() | |
# m_ProcessingHeader = ProcessingHeader() | |||
|
r487 | ||
systemHeaderObj = SystemHeader() | |||
r889 | |||
|
r487 | radarControllerHeaderObj = RadarControllerHeader() | |
# data = None | |||
r889 | |||
|
r487 | type = None | |
r889 | |||
|
r1092 | 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 | |||
|
r1092 | flagDecodeData = False # asumo q la data no esta decodificada | |
r889 | |||
|
r1092 | 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 | |||
|
r1092 | # Speed of ligth | |
|
r487 | 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 | |||
|
r1176 | error = (0, '') | |
def __str__(self): | |||
return '{} - {}'.format(self.type, self.getDatatime()) | |||
|
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 | |||
|
r1167 | return list(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: | |
|
r1092 | 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) | |
|
r1092 | datatime.append(self.ltctime + self.timeInterval + 1) | |
r889 | |||
|
r487 | datatime = numpy.array(datatime) | |
r889 | |||
|
r487 | return datatime | |
r889 | |||
|
r765 | def getFmaxTimeResponse(self): | |
r889 | |||
|
r1092 | period = (10**-6) * self.getDeltaH() / (0.15) | |
r889 | |||
|
r1092 | PRF = 1. / (period * self.nCohInt) | |
r889 | |||
|
r765 | fmax = PRF | |
r889 | |||
|
r765 | return fmax | |
r889 | |||
|
r487 | def getFmax(self): | |
|
r1092 | PRF = 1. / (self.ippSeconds * self.nCohInt) | |
r889 | |||
|
r765 | fmax = PRF | |
|
r487 | return fmax | |
r889 | |||
|
r487 | def getVmax(self): | |
r889 | |||
|
r1092 | _lambda = self.C / self.frequency | |
r889 | |||
|
r1092 | 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.") | |
|
r1092 | channelIndexList = property( | |
getChannelIndexList, "I'm the 'channelIndexList' property.") | |||
|
r487 | 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 | |||
|
r1092 | ||
|
r487 | class Voltage(JROData): | |
r889 | |||
|
r1092 | # data es un numpy array de 2 dmensiones (canales, alturas) | |
|
r487 | 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 | |||
|
r1092 | self.flagDecodeData = False # asumo q la data no esta decodificada | |
r889 | |||
|
r1092 | self.flagDeflipData = False # asumo q la data no esta sin flip | |
r889 | |||
|
r487 | self.flagShiftFFT = False | |
r889 | |||
|
r1092 | self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil | |
r889 | |||
|
r534 | self.profileIndex = 0 | |
r889 | |||
|
r1092 | 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: | |||
|
r1092 | daux = power[thisChannel, :].real | |
|
r568 | noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt) | |
r889 | |||
|
r568 | return noise | |
r889 | |||
|
r1092 | def getNoise(self, type=1, channel=None): | |
r889 | |||
|
r487 | if type == 1: | |
|
r568 | noise = self.getNoisebyHildebrand(channel) | |
r889 | |||
|
r720 | return noise | |
r889 | |||
|
r1092 | 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) | |
|
r1092 | powerdB = 10 * numpy.log10(power.real) | |
|
r850 | 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 | |||
|
r1092 | ||
|
r487 | class Spectra(JROData): | |
r889 | |||
|
r1092 | # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas) | |
|
r487 | data_spc = None | |
r889 | |||
|
r1092 | # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas) | |
|
r487 | data_cspc = None | |
r889 | |||
|
r1092 | # data dc es un numpy array de 2 dmensiones (canales, alturas) | |
|
r487 | data_dc = None | |
r889 | |||
|
r1092 | # data power | |
|
r832 | data_pwr = None | |
r889 | |||
|
r487 | nFFTPoints = None | |
r889 | |||
|
r487 | # nPairs = None | |
r889 | |||
|
r487 | pairsList = None | |
r889 | |||
|
r487 | nIncohInt = None | |
r889 | |||
|
r1092 | wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia | |
r889 | |||
|
r1092 | 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 | |||
|
r1092 | self.flagDecodeData = False # asumo q la data no esta decodificada | |
r889 | |||
|
r1092 | 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): | |
|
r1092 | 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: | |
|
r1092 | # this was estimated by getNoise Operation defined in jroproc_spectra.py | |
return self.noise_estimation | |||
|
r487 | else: | |
|
r1092 | noise = self.getNoisebyHildebrand( | |
xmin_index, xmax_index, ymin_index, ymax_index) | |||
|
r487 | return noise | |
r889 | |||
|
r765 | def getFreqRangeTimeResponse(self, extrapoints=0): | |
r889 | |||
|
r1092 | 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 | |||
|
r1092 | 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 | |||
|
r1092 | 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): | |||
r889 | |||
|
r1092 | deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor) | |
velrange = deltav * (numpy.arange(self.nFFTPoints + | |||
extrapoints) - self.nFFTPoints / 2.) # - deltav/2 | |||
r889 | |||
|
r487 | return velrange | |
r889 | |||
|
r487 | def getNPairs(self): | |
r889 | |||
|
r487 | return len(self.pairsList) | |
r889 | |||
|
r487 | def getPairsIndexList(self): | |
r889 | |||
|
r1167 | return list(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 | |||
|
r1092 | 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 | |||
|
r1160 | timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor | |
r889 | |||
|
r526 | return timeInterval | |
r889 | |||
|
r832 | def getPower(self): | |
r889 | |||
|
r832 | factor = self.normFactor | |
|
r1092 | z = self.data_spc / factor | |
r889 | z = numpy.where(numpy.isfinite(z), z, numpy.NAN) | ||
|
r832 | avg = numpy.average(z, axis=1) | |
r889 | |||
|
r1092 | 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: | |||
|
r1167 | raise ValueError("Pair %s is not in dataOut.pairsList" % ( | |
pair)) | |||
|
r1092 | pairsIndexList.append(self.pairsList.index(pair)) | |
r889 | for i in range(len(pairsIndexList)): | ||
pair = self.pairsList[pairsIndexList[i]] | |||
|
r1092 | ccf = numpy.average( | |
self.data_cspc[pairsIndexList[i], :, :], axis=0) | |||
r889 | powa = numpy.average(self.data_spc[pair[0], :, :], axis=0) | ||
powb = numpy.average(self.data_spc[pair[1], :, :], axis=0) | |||
|
r1092 | avgcoherenceComplex = ccf / numpy.sqrt(powa * powb) | |
r889 | if phase: | ||
data = numpy.arctan2(avgcoherenceComplex.imag, | |||
|
r1092 | avgcoherenceComplex.real) * 180 / numpy.pi | |
r889 | else: | ||
data = numpy.abs(avgcoherenceComplex) | |||
z.append(data) | |||
return numpy.array(z) | |||
|
r624 | def setValue(self, value): | |
r889 | |||
|
r1167 | print("This property should not be initialized") | |
r889 | |||
|
r624 | return | |
r889 | |||
|
r624 | nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.") | |
|
r1092 | pairsIndexList = property( | |
getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.") | |||
normFactor = property(getNormFactor, setValue, | |||
"I'm the 'getNormFactor' property.") | |||
|
r624 | flag_cspc = property(getFlagCspc, setValue) | |
flag_dc = property(getFlagDc, setValue) | |||
noise = property(getNoise, setValue, "I'm the 'nHeights' property.") | |||
|
r1092 | 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 | |||
|
r1092 | 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 | ||
|
r1092 | ||
|
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 | |||
|
r1092 | # Speed of ligth | |
|
r487 | C = 3e8 | |
r889 | |||
|
r487 | frequency = 49.92e6 | |
r889 | |||
|
r487 | realtime = False | |
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: | |
|
r1092 | 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 | |||
|
r1167 | return list(range(self.nChannels)) | |
r889 | |||
|
r1092 | 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.") | |||
|
r1092 | channelIndexList = property( | |
getChannelIndexList, "I'm the 'channelIndexList' property.") | |||
|
r487 | 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 | |
|
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 | |||
|
r1092 | self.flagDecodeData = False # asumo q la data no esta decodificada | |
r889 | |||
|
r1092 | 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 | |||
|
r1092 | def getNoise(self, mode=2): | |
r889 | |||
|
r502 | indR = numpy.where(self.lagR == 0)[0][0] | |
indT = numpy.where(self.lagT == 0)[0][0] | |||
r889 | |||
|
r1092 | jspectra0 = self.data_corr[:, :, indR, :] | |
|
r502 | jspectra = copy.copy(jspectra0) | |
r889 | |||
|
r502 | num_chan = jspectra.shape[0] | |
num_hei = jspectra.shape[2] | |||
r889 | |||
|
r1092 | freq_dc = jspectra.shape[1] / 2 | |
ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc | |||
r889 | |||
|
r1092 | if ind_vel[0] < 0: | |
|
r1167 | ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof | |
r889 | |||
if mode == 1: | |||
|
r1092 | jspectra[:, freq_dc, :] = ( | |
jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION | |||
r889 | |||
|
r502 | if mode == 2: | |
r889 | |||
|
r1092 | vel = numpy.array([-2, -1, 1, 2]) | |
xx = numpy.zeros([4, 4]) | |||
r889 | |||
|
r502 | for fil in range(4): | |
|
r1167 | xx[fil, :] = vel[fil]**numpy.asarray(list(range(4))) | |
r889 | |||
|
r502 | xx_inv = numpy.linalg.inv(xx) | |
|
r1092 | xx_aux = xx_inv[0, :] | |
r889 | |||
|
r502 | for ich in range(num_chan): | |
|
r1092 | yy = jspectra[ich, ind_vel, :] | |
jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy) | |||
|
r502 | ||
|
r1092 | junkid = jspectra[ich, freq_dc, :] <= 0 | |
|
r502 | cjunkid = sum(junkid) | |
r889 | |||
|
r502 | if cjunkid.any(): | |
|
r1092 | jspectra[ich, freq_dc, junkid.nonzero()] = ( | |
jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2 | |||
r889 | |||
|
r1092 | 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 | |||
|
r1092 | # 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) | |||
|
r1092 | 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 | |||
|
r1092 | 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 | |||
|
r1092 | ||
r922 | class Parameters(Spectra): | ||
|
r502 | ||
|
r1092 | experimentInfo = None # Information about the experiment | |
r889 | |||
|
r1092 | # Information from previous data | |
r889 | |||
|
r1092 | inputUnit = None # Type of data to be processed | |
r889 | |||
|
r1092 | operation = None # Type of operation to parametrize | |
r889 | |||
|
r1092 | # normFactor = None #Normalization Factor | |
r889 | |||
|
r1092 | groupList = None # List of Pairs, Groups, etc | |
r889 | |||
|
r1092 | # Parameters | |
r889 | |||
|
r1092 | data_param = None # Parameters obtained | |
r889 | |||
|
r1092 | data_pre = None # Data Pre Parametrization | |
r889 | |||
|
r1092 | data_SNR = None # Signal to Noise Ratio | |
r889 | |||
|
r608 | # heightRange = None #Heights | |
r889 | |||
|
r1092 | abscissaList = None # Abscissa, can be velocities, lags or time | |
r889 | |||
r923 | # noise = None #Noise Potency | ||
r889 | |||
|
r1092 | utctimeInit = None # Initial UTC time | |
r889 | |||
|
r1092 | paramInterval = None # Time interval to calculate Parameters in seconds | |
r889 | |||
|
r802 | useLocalTime = True | |
r889 | |||
|
r1092 | # Fitting | |
r889 | |||
|
r1092 | data_error = None # Error of the estimation | |
r889 | |||
constants = None | |||
|
r513 | library = None | |
r889 | |||
|
r1092 | # Output signal | |
r889 | |||
|
r1092 | outputInterval = None # Time interval to calculate output signal in seconds | |
r889 | |||
|
r1092 | data_output = None # Out signal | |
r889 | |||
|
r850 | nAvg = None | |
r889 | |||
r922 | noise_estimation = None | ||
|
r1092 | GauSPC = None # Fit gaussian SPC | |
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: | |
|
r1092 | time1 = self.utctimeInit - self.timeZone * 60 | |
|
r608 | else: | |
|
r802 | time1 = self.utctimeInit | |
r889 | |||
|
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): | |
|
r1167 | print("This property should not be initialized") | |
|
r1001 | ||
return | |||
|
r928 | def getNoise(self): | |
return self.spc_noise | |||
timeInterval = property(getTimeInterval) | |||
|
r1167 | noise = property(getNoise, setValue, "I'm the 'Noise' property.") |