jrodata.py
689 lines
| 16.5 KiB
| text/x-python
|
PythonLexer
|
r174 | ''' | |
$Author: murco $ | |||
$Id: JROData.py 173 2012-11-20 15:06:21Z murco $ | |||
''' | |||
import os, sys | |||
import copy | |||
import numpy | |||
|
r224 | import datetime | |
|
r174 | ||
|
r178 | from jroheaderIO import SystemHeader, RadarControllerHeader | |
|
r174 | ||
|
r199 | ||
|
r450 | def hildebrand_sekhon(data, navg): | |
|
r201 | ||
data = data.copy() | |||
|
r450 | sortdata = numpy.sort(data,axis=None) | |
lenOfData = len(sortdata) | |||
|
r199 | nums_min = lenOfData/10 | |
|
r450 | if (lenOfData/10) > 2: | |
|
r199 | nums_min = lenOfData/10 | |
else: | |||
|
r450 | nums_min = 2 | |
|
r445 | sump = 0. | |
|
r199 | ||
sumq = 0. | |||
j = 0 | |||
cont = 1 | |||
while((cont==1)and(j<lenOfData)): | |||
|
r445 | sump += sortdata[j] | |
|
r199 | ||
sumq += sortdata[j]**2 | |||
j += 1 | |||
if j > nums_min: | |||
|
r450 | rtest = float(j)/(j-1) + 1.0/navg | |
if ((sumq*j) > (rtest*sump**2)): | |||
|
r199 | j = j - 1 | |
|
r445 | sump = sump - sortdata[j] | |
sumq = sumq - sortdata[j]**2 | |||
|
r199 | cont = 0 | |
|
r450 | ||
lnoise = sump /j | |||
stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1)) | |||
|
r199 | return lnoise | |
|
r174 | class JROData: | |
# m_BasicHeader = BasicHeader() | |||
# m_ProcessingHeader = ProcessingHeader() | |||
systemHeaderObj = SystemHeader() | |||
radarControllerHeaderObj = RadarControllerHeader() | |||
# data = None | |||
type = None | |||
dtype = None | |||
|
r200 | # nChannels = None | |
|
r174 | ||
|
r207 | # nHeights = None | |
|
r174 | ||
nProfiles = None | |||
heightList = None | |||
channelList = None | |||
flagNoData = True | |||
flagTimeBlock = False | |||
|
r344 | useLocalTime = False | |
|
r174 | utctime = None | |
|
r344 | timeZone = None | |
dstFlag = None | |||
errorCount = None | |||
|
r174 | blocksize = None | |
nCode = None | |||
nBaud = None | |||
code = None | |||
|
r232 | flagDecodeData = False #asumo q la data no esta decodificada | |
|
r174 | ||
|
r232 | flagDeflipData = False #asumo q la data no esta sin flip | |
|
r174 | ||
flagShiftFFT = False | |||
ippSeconds = None | |||
timeInterval = None | |||
|
r199 | ||
nCohInt = None | |||
noise = None | |||
|
r201 | ||
|
r245 | windowOfFilter = 1 | |
|
r201 | #Speed of ligth | |
C = 3e8 | |||
frequency = 49.92e6 | |||
|
r360 | ||
realtime = False | |||
|
r453 | ||
beacon_heiIndexList = None | |||
|
r457 | ||
last_block = None | |||
blocknow = None | |||
|
r174 | ||
def __init__(self): | |||
raise ValueError, "This class has not been implemented" | |||
def copy(self, inputObj=None): | |||
if inputObj == None: | |||
return copy.deepcopy(self) | |||
for key in inputObj.__dict__.keys(): | |||
self.__dict__[key] = inputObj.__dict__[key] | |||
def deepcopy(self): | |||
return copy.deepcopy(self) | |||
|
r185 | def isEmpty(self): | |
|
r189 | return self.flagNoData | |
|
r185 | ||
|
r207 | def getNoise(self): | |
raise ValueError, "Not implemented" | |||
|
r200 | def getNChannels(self): | |
return len(self.channelList) | |||
def getChannelIndexList(self): | |||
return range(self.nChannels) | |||
|
r207 | def getNHeights(self): | |
return len(self.heightList) | |||
def getHeiRange(self, extrapoints=0): | |||
heis = self.heightList | |||
# deltah = self.heightList[1] - self.heightList[0] | |||
# | |||
# heis.append(self.heightList[-1]) | |||
return heis | |||
|
r344 | def getltctime(self): | |
if self.useLocalTime: | |||
return self.utctime - self.timeZone*60 | |||
return self.utctime | |||
|
r207 | def getDatatime(self): | |
|
r344 | datatime = datetime.datetime.utcfromtimestamp(self.ltctime) | |
|
r224 | return datatime | |
def getTimeRange(self): | |||
|
r207 | datatime = [] | |
|
r344 | datatime.append(self.ltctime) | |
datatime.append(self.ltctime + self.timeInterval) | |||
|
r207 | ||
datatime = numpy.array(datatime) | |||
return datatime | |||
def getFmax(self): | |||
PRF = 1./(self.ippSeconds * self.nCohInt) | |||
fmax = PRF/2. | |||
return fmax | |||
def getVmax(self): | |||
_lambda = self.C/self.frequency | |||
|
r227 | vmax = self.getFmax() * _lambda | |
|
r207 | ||
return vmax | |||
|
r200 | nChannels = property(getNChannels, "I'm the 'nChannel' property.") | |
channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.") | |||
|
r207 | nHeights = property(getNHeights, "I'm the 'nHeights' property.") | |
noise = property(getNoise, "I'm the 'nHeights' property.") | |||
|
r224 | datatime = property(getDatatime, "I'm the 'datatime' property") | |
|
r344 | ltctime = property(getltctime, "I'm the 'ltctime' property") | |
|
r207 | ||
|
r174 | class Voltage(JROData): | |
#data es un numpy array de 2 dmensiones (canales, alturas) | |||
data = None | |||
def __init__(self): | |||
''' | |||
Constructor | |||
''' | |||
self.radarControllerHeaderObj = RadarControllerHeader() | |||
self.systemHeaderObj = SystemHeader() | |||
self.type = "Voltage" | |||
self.data = None | |||
self.dtype = None | |||
|
r200 | # self.nChannels = 0 | |
|
r174 | ||
|
r207 | # self.nHeights = 0 | |
|
r174 | ||
self.nProfiles = None | |||
self.heightList = None | |||
self.channelList = None | |||
|
r200 | # self.channelIndexList = None | |
|
r174 | ||
self.flagNoData = True | |||
self.flagTimeBlock = False | |||
self.utctime = None | |||
|
r344 | self.timeZone = None | |
self.dstFlag = None | |||
self.errorCount = None | |||
|
r174 | self.nCohInt = None | |
self.blocksize = None | |||
|
r232 | ||
self.flagDecodeData = False #asumo q la data no esta decodificada | |||
self.flagDeflipData = False #asumo q la data no esta sin flip | |||
self.flagShiftFFT = False | |||
|
r199 | ||
def getNoisebyHildebrand(self): | |||
""" | |||
Determino el nivel de ruido usando el metodo Hildebrand-Sekhon | |||
Return: | |||
noiselevel | |||
""" | |||
|
r174 | ||
|
r199 | for channel in range(self.nChannels): | |
daux = self.data_spc[channel,:,:] | |||
self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt) | |||
return self.noise | |||
def getNoise(self, type = 1): | |||
self.noise = numpy.zeros(self.nChannels) | |||
if type == 1: | |||
noise = self.getNoisebyHildebrand() | |||
return 10*numpy.log10(noise) | |||
|
r174 | class Spectra(JROData): | |
#data es un numpy array de 2 dmensiones (canales, perfiles, alturas) | |||
data_spc = None | |||
#data es un numpy array de 2 dmensiones (canales, pares, alturas) | |||
data_cspc = None | |||
#data es un numpy array de 2 dmensiones (canales, alturas) | |||
data_dc = None | |||
nFFTPoints = None | |||
nPairs = None | |||
pairsList = None | |||
nIncohInt = None | |||
wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia | |||
nCohInt = None #se requiere para determinar el valor de timeInterval | |||
|
r445 | ippFactor = None | |
|
r174 | def __init__(self): | |
''' | |||
Constructor | |||
''' | |||
self.radarControllerHeaderObj = RadarControllerHeader() | |||
self.systemHeaderObj = SystemHeader() | |||
self.type = "Spectra" | |||
# self.data = None | |||
self.dtype = None | |||
|
r200 | # self.nChannels = 0 | |
|
r174 | ||
|
r207 | # self.nHeights = 0 | |
|
r174 | ||
self.nProfiles = None | |||
self.heightList = None | |||
self.channelList = None | |||
|
r200 | # self.channelIndexList = None | |
|
r174 | ||
self.flagNoData = True | |||
self.flagTimeBlock = False | |||
self.utctime = None | |||
|
r199 | self.nCohInt = None | |
|
r174 | self.nIncohInt = None | |
self.blocksize = None | |||
self.nFFTPoints = None | |||
self.wavelength = None | |||
|
r232 | ||
self.flagDecodeData = False #asumo q la data no esta decodificada | |||
self.flagDeflipData = False #asumo q la data no esta sin flip | |||
self.flagShiftFFT = False | |||
|
r448 | ||
self.ippFactor = 1 | |||
|
r450 | ||
self.noise = None | |||
|
r453 | ||
self.beacon_heiIndexList = [] | |||
|
r199 | ||
def getNoisebyHildebrand(self): | |||
""" | |||
Determino el nivel de ruido usando el metodo Hildebrand-Sekhon | |||
Return: | |||
noiselevel | |||
""" | |||
|
r452 | noise = numpy.zeros(self.nChannels) | |
|
r199 | for channel in range(self.nChannels): | |
daux = self.data_spc[channel,:,:] | |||
|
r452 | noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) | |
|
r199 | ||
|
r452 | return noise | |
|
r199 | ||
|
r451 | def getNoise(self): | |
|
r452 | if self.noise != None: | |
return self.noise | |||
else: | |||
noise = self.getNoisebyHildebrand() | |||
return noise | |||
|
r215 | ||
def getFreqRange(self, extrapoints=0): | |||
|
r445 | deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor) | |
|
r227 | freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2 | |
|
r215 | ||
return freqrange | |||
def getVelRange(self, extrapoints=0): | |||
|
r445 | deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor) | |
|
r215 | velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2 | |
return velrange | |||
def getNPairs(self): | |||
return len(self.pairsList) | |||
def getPairsIndexList(self): | |||
return range(self.nPairs) | |||
|
r245 | def getNormFactor(self): | |
pwcode = 1 | |||
if self.flagDecodeData: | |||
pwcode = numpy.sum(self.code[0]**2) | |||
|
r463 | #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter | |
normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter | |||
|
r245 | ||
return normFactor | |||
|
r266 | def getFlagCspc(self): | |
if self.data_cspc == None: | |||
return True | |||
return False | |||
def getFlagDc(self): | |||
if self.data_dc == None: | |||
return True | |||
return False | |||
|
r215 | nPairs = property(getNPairs, "I'm the 'nPairs' property.") | |
pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.") | |||
|
r245 | normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.") | |
|
r266 | flag_cspc = property(getFlagCspc) | |
flag_dc = property(getFlagDc) | |||
|
r174 | ||
class SpectraHeis(JROData): | |||
data_spc = None | |||
data_cspc = None | |||
data_dc = None | |||
nFFTPoints = None | |||
nPairs = None | |||
pairsList = None | |||
nIncohInt = None | |||
def __init__(self): | |||
self.radarControllerHeaderObj = RadarControllerHeader() | |||
self.systemHeaderObj = SystemHeader() | |||
self.type = "SpectraHeis" | |||
self.dtype = None | |||
|
r200 | # self.nChannels = 0 | |
|
r174 | ||
|
r207 | # self.nHeights = 0 | |
|
r174 | ||
self.nProfiles = None | |||
self.heightList = None | |||
self.channelList = None | |||
|
r200 | # self.channelIndexList = None | |
|
r174 | ||
self.flagNoData = True | |||
self.flagTimeBlock = False | |||
self.nPairs = 0 | |||
self.utctime = None | |||
self.blocksize = None | |||
|
r353 | ||
class Fits: | |||
|
r439 | heightList = None | |
channelList = None | |||
flagNoData = True | |||
flagTimeBlock = False | |||
useLocalTime = False | |||
utctime = None | |||
timeZone = None | |||
ippSeconds = None | |||
timeInterval = None | |||
nCohInt = None | |||
nIncohInt = None | |||
noise = None | |||
windowOfFilter = 1 | |||
#Speed of ligth | |||
C = 3e8 | |||
frequency = 49.92e6 | |||
realtime = False | |||
|
r353 | def __init__(self): | |
|
r439 | ||
self.type = "Fits" | |||
self.nProfiles = None | |||
self.heightList = None | |||
self.channelList = None | |||
# self.channelIndexList = None | |||
self.flagNoData = True | |||
|
r353 | self.utctime = None | |
|
r439 | ||
self.nCohInt = None | |||
self.nIncohInt = None | |||
self.useLocalTime = True | |||
# 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 = '' | |||
# | |||
|
r353 | ||
def getltctime(self): | |||
if self.useLocalTime: | |||
return self.utctime - self.timeZone*60 | |||
return self.utctime | |||
def getDatatime(self): | |||
datatime = datetime.datetime.utcfromtimestamp(self.ltctime) | |||
return datatime | |||
def getTimeRange(self): | |||
datatime = [] | |||
datatime.append(self.ltctime) | |||
datatime.append(self.ltctime + self.timeInterval) | |||
datatime = numpy.array(datatime) | |||
return datatime | |||
|
r439 | def getHeiRange(self): | |
heis = self.heightList | |||
return heis | |||
|
r353 | def isEmpty(self): | |
return self.flagNoData | |||
|
r439 | def getNHeights(self): | |
return len(self.heightList) | |||
def getNChannels(self): | |||
return len(self.channelList) | |||
def getChannelIndexList(self): | |||
return range(self.nChannels) | |||
def getNoise(self, type = 1): | |||
self.noise = numpy.zeros(self.nChannels) | |||
if type == 1: | |||
noise = self.getNoisebyHildebrand() | |||
if type == 2: | |||
noise = self.getNoisebySort() | |||
if type == 3: | |||
noise = self.getNoisebyWindow() | |||
return noise | |||
|
r353 | datatime = property(getDatatime, "I'm the 'datatime' property") | |
|
r439 | 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.") | |||
datatime = property(getDatatime, "I'm the 'datatime' property") | |||
ltctime = property(getltctime, "I'm the 'ltctime' property") | |||
|
r473 | ltctime = property(getltctime, "I'm the 'ltctime' property") | |
class AMISR: | |||
def __init__(self): | |||
self.flagNoData = True | |||
self.data = None | |||
self.utctime = None | |||
self.type = "AMISR" | |||
|
r474 | #propiedades para compatibilidad con Voltages | |
self.timeZone = 0#self.dataIn.timeZone | |||
self.dstFlag = 0#self.dataIn.dstFlag | |||
self.errorCount = 0#self.dataIn.errorCount | |||
self.useLocalTime = True#self.dataIn.useLocalTime | |||
self.radarControllerHeaderObj = None#self.dataIn.radarControllerHeaderObj.copy() | |||
self.systemHeaderObj = None#self.dataIn.systemHeaderObj.copy() | |||
self.channelList = [1]#self.dataIn.channelList esto solo aplica para el caso de AMISR | |||
self.dtype = numpy.dtype([('real','<f4'),('imag','<f4')]) | |||
self.flagTimeBlock = None#self.dataIn.flagTimeBlock | |||
#self.utctime = #self.firstdatatime | |||
self.flagDecodeData = None#self.dataIn.flagDecodeData #asumo q la data esta decodificada | |||
self.flagDeflipData = None#self.dataIn.flagDeflipData #asumo q la data esta sin flip | |||
self.nCohInt = 1#self.dataIn.nCohInt | |||
self.nIncohInt = 1 | |||
self.ippSeconds = 0.004#self.dataIn.ippSeconds, segun el filename/Setup/Tufile | |||
self.windowOfFilter = None#self.dataIn.windowOfFilter | |||
self.timeInterval = None#self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt | |||
self.frequency = 20000#self.dataIn.frequency | |||
self.realtime = 0#self.dataIn.realtime | |||
#actualizar en la lectura de datos | |||
self.heightList = None#self.dataIn.heightList | |||
self.nProfiles = None#self.dataOut.nFFTPoints | |||
self.nBaud = None#self.dataIn.nBaud | |||
self.nCode = None#self.dataIn.nCode | |||
self.code = None#self.dataIn.code | |||
|
r473 | def isEmpty(self): | |
return self.flagNoData |