##// END OF EJS Templates
schain
RSS
para pruebas em PC AMISR
para pruebas em PC AMISR
joabAM - - Load All Authors

File last commit:

r1541:77e2d8d2b2ef
r1541:77e2d8d2b2ef
Show More
jroproc_voltage.py
2361 lines | 77.1 KiB | text/x-python | PythonLexer
/ schainpy / model / proc / jroproc_voltage.py
History | Annotation | Raw |Copy content |Copy permalink
import sys
import numpy,math
from scipy import interpolate
from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
from schainpy.utils import log
from schainpy.model.io.utils import getHei_index
from time import time
#import datetime
import numpy
#import copy
from schainpy.model.data import _noise
class VoltageProc(ProcessingUnit):
def __init__(self):
ProcessingUnit.__init__(self)
self.dataOut = Voltage()
self.flip = 1
self.setupReq = False
def run(self):
#print("running volt proc")
if self.dataIn.type == 'AMISR':
self.__updateObjFromAmisrInput()
if self.dataOut.buffer_empty:
if self.dataIn.type == 'Voltage':
self.dataOut.copy(self.dataIn)
#print("new volts reading")
def __updateObjFromAmisrInput(self):
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.flagNoData = self.dataIn.flagNoData
self.dataOut.data = self.dataIn.data
self.dataOut.utctime = self.dataIn.utctime
self.dataOut.channelList = self.dataIn.channelList
#self.dataOut.timeInterval = self.dataIn.timeInterval
self.dataOut.heightList = self.dataIn.heightList
self.dataOut.nProfiles = self.dataIn.nProfiles
self.dataOut.nCohInt = self.dataIn.nCohInt
self.dataOut.ippSeconds = self.dataIn.ippSeconds
self.dataOut.frequency = self.dataIn.frequency
self.dataOut.azimuth = self.dataIn.azimuth
self.dataOut.zenith = self.dataIn.zenith
self.dataOut.beam.codeList = self.dataIn.beam.codeList
self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
class selectChannels(Operation):
def run(self, dataOut, channelList=None):
self.channelList = channelList
if self.channelList == None:
print("Missing channelList")
return dataOut
channelIndexList = []
if type(dataOut.channelList) is not list: #leer array desde HDF5
try:
dataOut.channelList = dataOut.channelList.tolist()
except Exception as e:
print("Select Channels: ",e)
for channel in self.channelList:
if channel not in dataOut.channelList:
raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
index = dataOut.channelList.index(channel)
channelIndexList.append(index)
dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
return dataOut
def selectChannelsByIndex(self, dataOut, channelIndexList):
"""
Selecciona un bloque de datos en base a canales segun el channelIndexList
Input:
channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
Affected:
dataOut.data
dataOut.channelIndexList
dataOut.nChannels
dataOut.m_ProcessingHeader.totalSpectra
dataOut.systemHeaderObj.numChannels
dataOut.m_ProcessingHeader.blockSize
Return:
None
"""
#print("selectChannelsByIndex")
# for channelIndex in channelIndexList:
# if channelIndex not in dataOut.channelIndexList:
# raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
if dataOut.type == 'Voltage':
if dataOut.flagDataAsBlock:
"""
Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
"""
data = dataOut.data[channelIndexList,:,:]
else:
data = dataOut.data[channelIndexList,:]
dataOut.data = data
# dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
dataOut.channelList = range(len(channelIndexList))
elif dataOut.type == 'Spectra':
if hasattr(dataOut, 'data_spc'):
if dataOut.data_spc is None:
raise ValueError("data_spc is None")
return dataOut
else:
data_spc = dataOut.data_spc[channelIndexList, :]
dataOut.data_spc = data_spc
# if hasattr(dataOut, 'data_dc') :# and
# if dataOut.data_dc is None:
# raise ValueError("data_dc is None")
# return dataOut
# else:
# data_dc = dataOut.data_dc[channelIndexList, :]
# dataOut.data_dc = data_dc
# dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
dataOut.channelList = channelIndexList
dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
return dataOut
def __selectPairsByChannel(self, dataOut, channelList=None):
#print("__selectPairsByChannel")
if channelList == None:
return
pairsIndexListSelected = []
for pairIndex in dataOut.pairsIndexList:
# First pair
if dataOut.pairsList[pairIndex][0] not in channelList:
continue
# Second pair
if dataOut.pairsList[pairIndex][1] not in channelList:
continue
pairsIndexListSelected.append(pairIndex)
if not pairsIndexListSelected:
dataOut.data_cspc = None
dataOut.pairsList = []
return
dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
dataOut.pairsList = [dataOut.pairsList[i]
for i in pairsIndexListSelected]
return dataOut
class selectHeights(Operation):
def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
"""
Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
minHei <= height <= maxHei
Input:
minHei : valor minimo de altura a considerar
maxHei : valor maximo de altura a considerar
Affected:
Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
Return:
1 si el metodo se ejecuto con exito caso contrario devuelve 0
"""
self.dataOut = dataOut
if minHei and maxHei:
if (minHei < dataOut.heightList[0]):
minHei = dataOut.heightList[0]
if (maxHei > dataOut.heightList[-1]):
maxHei = dataOut.heightList[-1]
minIndex = 0
maxIndex = 0
heights = dataOut.heightList
inda = numpy.where(heights >= minHei)
indb = numpy.where(heights <= maxHei)
try:
minIndex = inda[0][0]
except:
minIndex = 0
try:
maxIndex = indb[0][-1]
except:
maxIndex = len(heights)
self.selectHeightsByIndex(minIndex, maxIndex)
return dataOut
def selectHeightsByIndex(self, minIndex, maxIndex):
"""
Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
minIndex <= index <= maxIndex
Input:
minIndex : valor de indice minimo de altura a considerar
maxIndex : valor de indice maximo de altura a considerar
Affected:
self.dataOut.data
self.dataOut.heightList
Return:
1 si el metodo se ejecuto con exito caso contrario devuelve 0
"""
if self.dataOut.type == 'Voltage':
if (minIndex < 0) or (minIndex > maxIndex):
raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
if (maxIndex >= self.dataOut.nHeights):
maxIndex = self.dataOut.nHeights
#voltage
if self.dataOut.flagDataAsBlock:
"""
Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
"""
data = self.dataOut.data[:,:, minIndex:maxIndex]
else:
data = self.dataOut.data[:, minIndex:maxIndex]
# firstHeight = self.dataOut.heightList[minIndex]
self.dataOut.data = data
self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
if self.dataOut.nHeights <= 1:
raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
elif self.dataOut.type == 'Spectra':
if (minIndex < 0) or (minIndex > maxIndex):
raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
minIndex, maxIndex))
if (maxIndex >= self.dataOut.nHeights):
maxIndex = self.dataOut.nHeights - 1
# Spectra
data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
data_cspc = None
if self.dataOut.data_cspc is not None:
data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
data_dc = None
if self.dataOut.data_dc is not None:
data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
self.dataOut.data_spc = data_spc
self.dataOut.data_cspc = data_cspc
self.dataOut.data_dc = data_dc
self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
return 1
class filterByHeights(Operation):
def run(self, dataOut, window):
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
if window == None:
window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
newdelta = deltaHeight * window
r = dataOut.nHeights % window
newheights = (dataOut.nHeights-r)/window
if newheights <= 1:
raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
if dataOut.flagDataAsBlock:
"""
Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
"""
buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
buffer = numpy.sum(buffer,3)
else:
buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
buffer = numpy.sum(buffer,2)
dataOut.data = buffer
dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
dataOut.windowOfFilter = window
return dataOut
class setH0(Operation):
def run(self, dataOut, h0, deltaHeight = None):
if not deltaHeight:
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
nHeights = dataOut.nHeights
newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
dataOut.heightList = newHeiRange
return dataOut
class deFlip(Operation):
def run(self, dataOut, channelList = []):
data = dataOut.data.copy()
if dataOut.flagDataAsBlock:
flip = self.flip
profileList = list(range(dataOut.nProfiles))
if not channelList:
for thisProfile in profileList:
data[:,thisProfile,:] = data[:,thisProfile,:]*flip
flip *= -1.0
else:
for thisChannel in channelList:
if thisChannel not in dataOut.channelList:
continue
for thisProfile in profileList:
data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
flip *= -1.0
self.flip = flip
else:
if not channelList:
data[:,:] = data[:,:]*self.flip
else:
for thisChannel in channelList:
if thisChannel not in dataOut.channelList:
continue
data[thisChannel,:] = data[thisChannel,:]*self.flip
self.flip *= -1.
dataOut.data = data
return dataOut
class setAttribute(Operation):
'''
Set an arbitrary attribute(s) to dataOut
'''
def __init__(self):
Operation.__init__(self)
self._ready = False
def run(self, dataOut, **kwargs):
for key, value in kwargs.items():
setattr(dataOut, key, value)
return dataOut
@MPDecorator
class printAttribute(Operation):
'''
Print an arbitrary attribute of dataOut
'''
def __init__(self):
Operation.__init__(self)
def run(self, dataOut, attributes):
if isinstance(attributes, str):
attributes = [attributes]
for attr in attributes:
if hasattr(dataOut, attr):
log.log(getattr(dataOut, attr), attr)
class interpolateHeights(Operation):
def run(self, dataOut, topLim, botLim):
#69 al 72 para julia
#82-84 para meteoros
if len(numpy.shape(dataOut.data))==2:
sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
#dataOut.data[:,botLim:limSup+1] = sampInterp
dataOut.data[:,botLim:topLim+1] = sampInterp
else:
nHeights = dataOut.data.shape[2]
x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
f = interpolate.interp1d(x, y, axis = 2)
xnew = numpy.arange(botLim,topLim+1)
ynew = f(xnew)
dataOut.data[:,:,botLim:topLim+1] = ynew
return dataOut
class CohInt(Operation):
isConfig = False
__profIndex = 0
__byTime = False
__initime = None
__lastdatatime = None
__integrationtime = None
__buffer = None
__bufferStride = []
__dataReady = False
__profIndexStride = 0
__dataToPutStride = False
n = None
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
"""
Set the parameters of the integration class.
Inputs:
n : Number of coherent integrations
timeInterval : Time of integration. If the parameter "n" is selected this one does not work
overlapping :
"""
self.__initime = None
self.__lastdatatime = 0
self.__buffer = None
self.__dataReady = False
self.byblock = byblock
self.stride = stride
if n == None and timeInterval == None:
raise ValueError("n or timeInterval should be specified ...")
if n != None:
self.n = n
self.__byTime = False
else:
self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
self.n = 9999
self.__byTime = True
if overlapping:
self.__withOverlapping = True
self.__buffer = None
else:
self.__withOverlapping = False
self.__buffer = 0
self.__profIndex = 0
def putData(self, data):
"""
Add a profile to the __buffer and increase in one the __profileIndex
"""
if not self.__withOverlapping:
self.__buffer += data.copy()
self.__profIndex += 1
return
#Overlapping data
nChannels, nHeis = data.shape
data = numpy.reshape(data, (1, nChannels, nHeis))
#If the buffer is empty then it takes the data value
if self.__buffer is None:
self.__buffer = data
self.__profIndex += 1
return
#If the buffer length is lower than n then stakcing the data value
if self.__profIndex < self.n:
self.__buffer = numpy.vstack((self.__buffer, data))
self.__profIndex += 1
return
#If the buffer length is equal to n then replacing the last buffer value with the data value
self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
self.__buffer[self.n-1] = data
self.__profIndex = self.n
return
def pushData(self):
"""
Return the sum of the last profiles and the profiles used in the sum.
Affected:
self.__profileIndex
"""
if not self.__withOverlapping:
data = self.__buffer
n = self.__profIndex
self.__buffer = 0
self.__profIndex = 0
return data, n
#Integration with Overlapping
data = numpy.sum(self.__buffer, axis=0)
# print data
# raise
n = self.__profIndex
return data, n
def byProfiles(self, data):
self.__dataReady = False
avgdata = None
# n = None
# print data
# raise
self.putData(data)
if self.__profIndex == self.n:
avgdata, n = self.pushData()
self.__dataReady = True
return avgdata
def byTime(self, data, datatime):
self.__dataReady = False
avgdata = None
n = None
self.putData(data)
if (datatime - self.__initime) >= self.__integrationtime:
avgdata, n = self.pushData()
self.n = n
self.__dataReady = True
return avgdata
def integrateByStride(self, data, datatime):
# print data
if self.__profIndex == 0:
self.__buffer = [[data.copy(), datatime]]
else:
self.__buffer.append([data.copy(),datatime])
self.__profIndex += 1
self.__dataReady = False
if self.__profIndex == self.n * self.stride :
self.__dataToPutStride = True
self.__profIndexStride = 0
self.__profIndex = 0
self.__bufferStride = []
for i in range(self.stride):
current = self.__buffer[i::self.stride]
data = numpy.sum([t[0] for t in current], axis=0)
avgdatatime = numpy.average([t[1] for t in current])
# print data
self.__bufferStride.append((data, avgdatatime))
if self.__dataToPutStride:
self.__dataReady = True
self.__profIndexStride += 1
if self.__profIndexStride == self.stride:
self.__dataToPutStride = False
# print self.__bufferStride[self.__profIndexStride - 1]
# raise
return self.__bufferStride[self.__profIndexStride - 1]
return None, None
def integrate(self, data, datatime=None):
if self.__initime == None:
self.__initime = datatime
if self.__byTime:
avgdata = self.byTime(data, datatime)
else:
avgdata = self.byProfiles(data)
self.__lastdatatime = datatime
if avgdata is None:
return None, None
avgdatatime = self.__initime
deltatime = datatime - self.__lastdatatime
if not self.__withOverlapping:
self.__initime = datatime
else:
self.__initime += deltatime
return avgdata, avgdatatime
def integrateByBlock(self, dataOut):
times = int(dataOut.data.shape[1]/self.n)
avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
id_min = 0
id_max = self.n
for i in range(times):
junk = dataOut.data[:,id_min:id_max,:]
avgdata[:,i,:] = junk.sum(axis=1)
id_min += self.n
id_max += self.n
timeInterval = dataOut.ippSeconds*self.n
avgdatatime = (times - 1) * timeInterval + dataOut.utctime
self.__dataReady = True
return avgdata, avgdatatime
def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
if not self.isConfig:
self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
self.isConfig = True
if dataOut.flagDataAsBlock:
"""
Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
"""
avgdata, avgdatatime = self.integrateByBlock(dataOut)
dataOut.nProfiles /= self.n
else:
if stride is None:
avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
else:
avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
# dataOut.timeInterval *= n
dataOut.flagNoData = True
if self.__dataReady:
dataOut.data = avgdata
if not dataOut.flagCohInt:
dataOut.nCohInt *= self.n
dataOut.flagCohInt = True
dataOut.utctime = avgdatatime
# print avgdata, avgdatatime
# raise
# dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
dataOut.flagNoData = False
return dataOut
class Decoder(Operation):
isConfig = False
__profIndex = 0
code = None
nCode = None
nBaud = None
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.times = None
self.osamp = None
# self.__setValues = False
self.isConfig = False
self.setupReq = False
def setup(self, code, osamp, dataOut):
self.__profIndex = 0
self.code = code
self.nCode = len(code)
self.nBaud = len(code[0])
if (osamp != None) and (osamp >1):
self.osamp = osamp
self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
self.nBaud = self.nBaud*self.osamp
self.__nChannels = dataOut.nChannels
self.__nProfiles = dataOut.nProfiles
self.__nHeis = dataOut.nHeights
if self.__nHeis < self.nBaud:
raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
#Frequency
__codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
__codeBuffer[:,0:self.nBaud] = self.code
self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
if dataOut.flagDataAsBlock:
self.ndatadec = self.__nHeis #- self.nBaud + 1
self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
else:
#Time
self.ndatadec = self.__nHeis #- self.nBaud + 1
self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
def __convolutionInFreq(self, data):
fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
fft_data = numpy.fft.fft(data, axis=1)
conv = fft_data*fft_code
data = numpy.fft.ifft(conv,axis=1)
return data
def __convolutionInFreqOpt(self, data):
raise NotImplementedError
def __convolutionInTime(self, data):
code = self.code[self.__profIndex]
for i in range(self.__nChannels):
self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
return self.datadecTime
def __convolutionByBlockInTime(self, data):
repetitions = int(self.__nProfiles / self.nCode)
junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
junk = junk.flatten()
code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
profilesList = range(self.__nProfiles)
for i in range(self.__nChannels):
for j in profilesList:
self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
return self.datadecTime
def __convolutionByBlockInFreq(self, data):
raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
fft_data = numpy.fft.fft(data, axis=2)
conv = fft_data*fft_code
data = numpy.fft.ifft(conv,axis=2)
return data
def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
if dataOut.flagDecodeData:
print("This data is already decoded, recoding again ...")
if not self.isConfig:
if code is None:
if dataOut.code is None:
raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
code = dataOut.code
else:
code = numpy.array(code).reshape(nCode,nBaud)
self.setup(code, osamp, dataOut)
self.isConfig = True
if mode == 3:
sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
if times != None:
sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
if self.code is None:
print("Fail decoding: Code is not defined.")
return
self.__nProfiles = dataOut.nProfiles
datadec = None
if mode == 3:
mode = 0
if dataOut.flagDataAsBlock:
"""
Decoding when data have been read as block,
"""
if mode == 0:
datadec = self.__convolutionByBlockInTime(dataOut.data)
if mode == 1:
datadec = self.__convolutionByBlockInFreq(dataOut.data)
else:
"""
Decoding when data have been read profile by profile
"""
if mode == 0:
datadec = self.__convolutionInTime(dataOut.data)
if mode == 1:
datadec = self.__convolutionInFreq(dataOut.data)
if mode == 2:
datadec = self.__convolutionInFreqOpt(dataOut.data)
if datadec is None:
raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
dataOut.code = self.code
dataOut.nCode = self.nCode
dataOut.nBaud = self.nBaud
dataOut.data = datadec
dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
dataOut.flagDecodeData = True #asumo q la data esta decodificada
if self.__profIndex == self.nCode-1:
self.__profIndex = 0
return dataOut
self.__profIndex += 1
return dataOut
# dataOut.flagDeflipData = True #asumo q la data no esta sin flip
class ProfileConcat(Operation):
isConfig = False
buffer = None
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.profileIndex = 0
def reset(self):
self.buffer = numpy.zeros_like(self.buffer)
self.start_index = 0
self.times = 1
def setup(self, data, m, n=1):
self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
self.nHeights = data.shape[1]#.nHeights
self.start_index = 0
self.times = 1
def concat(self, data):
self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
self.start_index = self.start_index + self.nHeights
def run(self, dataOut, m):
dataOut.flagNoData = True
if not self.isConfig:
self.setup(dataOut.data, m, 1)
self.isConfig = True
if dataOut.flagDataAsBlock:
raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
else:
self.concat(dataOut.data)
self.times += 1
if self.times > m:
dataOut.data = self.buffer
self.reset()
dataOut.flagNoData = False
# se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
dataOut.ippSeconds *= m
return dataOut
class ProfileSelector(Operation):
profileIndex = None
# Tamanho total de los perfiles
nProfiles = None
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.profileIndex = 0
def incProfileIndex(self):
self.profileIndex += 1
if self.profileIndex >= self.nProfiles:
self.profileIndex = 0
def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
if profileIndex < minIndex:
return False
if profileIndex > maxIndex:
return False
return True
def isThisProfileInList(self, profileIndex, profileList):
if profileIndex not in profileList:
return False
return True
def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
"""
ProfileSelector:
Inputs:
profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
"""
if rangeList is not None:
if type(rangeList[0]) not in (tuple, list):
rangeList = [rangeList]
dataOut.flagNoData = True
if dataOut.flagDataAsBlock:
"""
data dimension = [nChannels, nProfiles, nHeis]
"""
if profileList != None:
dataOut.data = dataOut.data[:,profileList,:]
if profileRangeList != None:
minIndex = profileRangeList[0]
maxIndex = profileRangeList[1]
profileList = list(range(minIndex, maxIndex+1))
dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
if rangeList != None:
profileList = []
for thisRange in rangeList:
minIndex = thisRange[0]
maxIndex = thisRange[1]
profileList.extend(list(range(minIndex, maxIndex+1)))
dataOut.data = dataOut.data[:,profileList,:]
dataOut.nProfiles = len(profileList)
dataOut.profileIndex = dataOut.nProfiles - 1
dataOut.flagNoData = False
return dataOut
"""
data dimension = [nChannels, nHeis]
"""
if profileList != None:
if self.isThisProfileInList(dataOut.profileIndex, profileList):
self.nProfiles = len(profileList)
dataOut.nProfiles = self.nProfiles
dataOut.profileIndex = self.profileIndex
dataOut.flagNoData = False
self.incProfileIndex()
return dataOut
if profileRangeList != None:
minIndex = profileRangeList[0]
maxIndex = profileRangeList[1]
if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
self.nProfiles = maxIndex - minIndex + 1
dataOut.nProfiles = self.nProfiles
dataOut.profileIndex = self.profileIndex
dataOut.flagNoData = False
self.incProfileIndex()
return dataOut
if rangeList != None:
nProfiles = 0
for thisRange in rangeList:
minIndex = thisRange[0]
maxIndex = thisRange[1]
nProfiles += maxIndex - minIndex + 1
for thisRange in rangeList:
minIndex = thisRange[0]
maxIndex = thisRange[1]
if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
self.nProfiles = nProfiles
dataOut.nProfiles = self.nProfiles
dataOut.profileIndex = self.profileIndex
dataOut.flagNoData = False
self.incProfileIndex()
break
return dataOut
if beam != None: #beam is only for AMISR data
if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
dataOut.flagNoData = False
dataOut.profileIndex = self.profileIndex
self.incProfileIndex()
return dataOut
raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
class Reshaper(Operation):
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.__buffer = None
self.__nitems = 0
def __appendProfile(self, dataOut, nTxs):
if self.__buffer is None:
shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
ini = dataOut.nHeights * self.__nitems
end = ini + dataOut.nHeights
self.__buffer[:, ini:end] = dataOut.data
self.__nitems += 1
return int(self.__nitems*nTxs)
def __getBuffer(self):
if self.__nitems == int(1./self.__nTxs):
self.__nitems = 0
return self.__buffer.copy()
return None
def __checkInputs(self, dataOut, shape, nTxs):
if shape is None and nTxs is None:
raise ValueError("Reshaper: shape of factor should be defined")
if nTxs:
if nTxs < 0:
raise ValueError("nTxs should be greater than 0")
if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
return shape, nTxs
if len(shape) != 2 and len(shape) != 3:
raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
if len(shape) == 2:
shape_tuple = [dataOut.nChannels]
shape_tuple.extend(shape)
else:
shape_tuple = list(shape)
nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
return shape_tuple, nTxs
def run(self, dataOut, shape=None, nTxs=None):
shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
dataOut.flagNoData = True
profileIndex = None
if dataOut.flagDataAsBlock:
dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
dataOut.flagNoData = False
profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
else:
if self.__nTxs < 1:
self.__appendProfile(dataOut, self.__nTxs)
new_data = self.__getBuffer()
if new_data is not None:
dataOut.data = new_data
dataOut.flagNoData = False
profileIndex = dataOut.profileIndex*nTxs
else:
raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
dataOut.profileIndex = profileIndex
dataOut.ippSeconds /= self.__nTxs
return dataOut
class SplitProfiles(Operation):
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
def run(self, dataOut, n):
dataOut.flagNoData = True
profileIndex = None
if dataOut.flagDataAsBlock:
#nchannels, nprofiles, nsamples
shape = dataOut.data.shape
if shape[2] % n != 0:
raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
new_shape = shape[0], shape[1]*n, int(shape[2]/n)
dataOut.data = numpy.reshape(dataOut.data, new_shape)
dataOut.flagNoData = False
profileIndex = int(dataOut.nProfiles/n) - 1
else:
raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
dataOut.nProfiles = int(dataOut.nProfiles*n)
dataOut.profileIndex = profileIndex
dataOut.ippSeconds /= n
return dataOut
class CombineProfiles(Operation):
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.__remData = None
self.__profileIndex = 0
def run(self, dataOut, n):
dataOut.flagNoData = True
profileIndex = None
if dataOut.flagDataAsBlock:
#nchannels, nprofiles, nsamples
shape = dataOut.data.shape
new_shape = shape[0], shape[1]/n, shape[2]*n
if shape[1] % n != 0:
raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
dataOut.data = numpy.reshape(dataOut.data, new_shape)
dataOut.flagNoData = False
profileIndex = int(dataOut.nProfiles*n) - 1
else:
#nchannels, nsamples
if self.__remData is None:
newData = dataOut.data
else:
newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
self.__profileIndex += 1
if self.__profileIndex < n:
self.__remData = newData
#continue
return
self.__profileIndex = 0
self.__remData = None
dataOut.data = newData
dataOut.flagNoData = False
profileIndex = dataOut.profileIndex/n
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
dataOut.nProfiles = int(dataOut.nProfiles/n)
dataOut.profileIndex = profileIndex
dataOut.ippSeconds *= n
return dataOut
class PulsePairVoltage(Operation):
'''
Function PulsePair(Signal Power, Velocity)
The real component of Lag[0] provides Intensity Information
The imag component of Lag[1] Phase provides Velocity Information
Configuration Parameters:
nPRF = Number of Several PRF
theta = Degree Azimuth angel Boundaries
Input:
self.dataOut
lag[N]
Affected:
self.dataOut.spc
'''
isConfig = False
__profIndex = 0
__initime = None
__lastdatatime = None
__buffer = None
noise = None
__dataReady = False
n = None
__nch = 0
__nHeis = 0
removeDC = False
ipp = None
lambda_ = 0
def __init__(self,**kwargs):
Operation.__init__(self,**kwargs)
def setup(self, dataOut, n = None, removeDC=False):
'''
n= Numero de PRF's de entrada
'''
self.__initime = None
self.__lastdatatime = 0
self.__dataReady = False
self.__buffer = 0
self.__profIndex = 0
self.noise = None
self.__nch = dataOut.nChannels
self.__nHeis = dataOut.nHeights
self.removeDC = removeDC
self.lambda_ = 3.0e8/(9345.0e6)
self.ippSec = dataOut.ippSeconds
self.nCohInt = dataOut.nCohInt
if n == None:
raise ValueError("n should be specified.")
if n != None:
if n<2:
raise ValueError("n should be greater than 2")
self.n = n
self.__nProf = n
self.__buffer = numpy.zeros((dataOut.nChannels,
n,
dataOut.nHeights),
dtype='complex')
def putData(self,data):
'''
Add a profile to he __buffer and increase in one the __profiel Index
'''
self.__buffer[:,self.__profIndex,:]= data
self.__profIndex += 1
return
def pushData(self,dataOut):
'''
Return the PULSEPAIR and the profiles used in the operation
Affected : self.__profileIndex
'''
#----------------- Remove DC-----------------------------------
if self.removeDC==True:
mean = numpy.mean(self.__buffer,1)
tmp = mean.reshape(self.__nch,1,self.__nHeis)
dc= numpy.tile(tmp,[1,self.__nProf,1])
self.__buffer = self.__buffer - dc
#------------------Calculo de Potencia ------------------------
pair0 = self.__buffer*numpy.conj(self.__buffer)
pair0 = pair0.real
lag_0 = numpy.sum(pair0,1)
#------------------Calculo de Ruido x canal--------------------
self.noise = numpy.zeros(self.__nch)
for i in range(self.__nch):
daux = numpy.sort(pair0[i,:,:],axis= None)
self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
self.noise = self.noise.reshape(self.__nch,1)
self.noise = numpy.tile(self.noise,[1,self.__nHeis])
noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
#------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
#------------------ P= S+N ,P=lag_0/N ---------------------------------
#-------------------- Power --------------------------------------------------
data_power = lag_0/(self.n*self.nCohInt)
#------------------ Senal ---------------------------------------------------
data_intensity = pair0 - noise_buffer
data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
#data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
for i in range(self.__nch):
for j in range(self.__nHeis):
if data_intensity[i][j] < 0:
data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
#----------------- Calculo de Frecuencia y Velocidad doppler--------
pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
lag_1 = numpy.sum(pair1,1)
data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
data_velocity = (self.lambda_/2.0)*data_freq
#---------------- Potencia promedio estimada de la Senal-----------
lag_0 = lag_0/self.n
S = lag_0-self.noise
#---------------- Frecuencia Doppler promedio ---------------------
lag_1 = lag_1/(self.n-1)
R1 = numpy.abs(lag_1)
#---------------- Calculo del SNR----------------------------------
data_snrPP = S/self.noise
for i in range(self.__nch):
for j in range(self.__nHeis):
if data_snrPP[i][j] < 1.e-20:
data_snrPP[i][j] = 1.e-20
#----------------- Calculo del ancho espectral ----------------------
L = S/R1
L = numpy.where(L<0,1,L)
L = numpy.log(L)
tmp = numpy.sqrt(numpy.absolute(L))
data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
n = self.__profIndex
self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
self.__profIndex = 0
return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
def pulsePairbyProfiles(self,dataOut):
self.__dataReady = False
data_power = None
data_intensity = None
data_velocity = None
data_specwidth = None
data_snrPP = None
self.putData(data=dataOut.data)
if self.__profIndex == self.n:
data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
self.__dataReady = True
return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
def pulsePairOp(self, dataOut, datatime= None):
if self.__initime == None:
self.__initime = datatime
data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
self.__lastdatatime = datatime
if data_power is None:
return None, None, None,None,None,None
avgdatatime = self.__initime
deltatime = datatime - self.__lastdatatime
self.__initime = datatime
return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
if not self.isConfig:
self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
self.isConfig = True
data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
dataOut.flagNoData = True
if self.__dataReady:
dataOut.nCohInt *= self.n
dataOut.dataPP_POW = data_intensity # S
dataOut.dataPP_POWER = data_power # P
dataOut.dataPP_DOP = data_velocity
dataOut.dataPP_SNR = data_snrPP
dataOut.dataPP_WIDTH = data_specwidth
dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
dataOut.utctime = avgdatatime
dataOut.flagNoData = False
return dataOut
# import collections
# from scipy.stats import mode
#
# class Synchronize(Operation):
#
# isConfig = False
# __profIndex = 0
#
# def __init__(self, **kwargs):
#
# Operation.__init__(self, **kwargs)
# # self.isConfig = False
# self.__powBuffer = None
# self.__startIndex = 0
# self.__pulseFound = False
#
# def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
#
# #Read data
#
# powerdB = dataOut.getPower(channel = channel)
# noisedB = dataOut.getNoise(channel = channel)[0]
#
# self.__powBuffer.extend(powerdB.flatten())
#
# dataArray = numpy.array(self.__powBuffer)
#
# filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
#
# maxValue = numpy.nanmax(filteredPower)
#
# if maxValue < noisedB + 10:
# #No se encuentra ningun pulso de transmision
# return None
#
# maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
#
# if len(maxValuesIndex) < 2:
# #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
# return None
#
# phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
#
# #Seleccionar solo valores con un espaciamiento de nSamples
# pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
#
# if len(pulseIndex) < 2:
# #Solo se encontro un pulso de transmision con ancho mayor a 1
# return None
#
# spacing = pulseIndex[1:] - pulseIndex[:-1]
#
# #remover senales que se distancien menos de 10 unidades o muestras
# #(No deberian existir IPP menor a 10 unidades)
#
# realIndex = numpy.where(spacing > 10 )[0]
#
# if len(realIndex) < 2:
# #Solo se encontro un pulso de transmision con ancho mayor a 1
# return None
#
# #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
# realPulseIndex = pulseIndex[realIndex]
#
# period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
#
# print "IPP = %d samples" %period
#
# self.__newNSamples = dataOut.nHeights #int(period)
# self.__startIndex = int(realPulseIndex[0])
#
# return 1
#
#
# def setup(self, nSamples, nChannels, buffer_size = 4):
#
# self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
# maxlen = buffer_size*nSamples)
#
# bufferList = []
#
# for i in range(nChannels):
# bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
# maxlen = buffer_size*nSamples)
#
# bufferList.append(bufferByChannel)
#
# self.__nSamples = nSamples
# self.__nChannels = nChannels
# self.__bufferList = bufferList
#
# def run(self, dataOut, channel = 0):
#
# if not self.isConfig:
# nSamples = dataOut.nHeights
# nChannels = dataOut.nChannels
# self.setup(nSamples, nChannels)
# self.isConfig = True
#
# #Append new data to internal buffer
# for thisChannel in range(self.__nChannels):
# bufferByChannel = self.__bufferList[thisChannel]
# bufferByChannel.extend(dataOut.data[thisChannel])
#
# if self.__pulseFound:
# self.__startIndex -= self.__nSamples
#
# #Finding Tx Pulse
# if not self.__pulseFound:
# indexFound = self.__findTxPulse(dataOut, channel)
#
# if indexFound == None:
# dataOut.flagNoData = True
# return
#
# self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
# self.__pulseFound = True
# self.__startIndex = indexFound
#
# #If pulse was found ...
# for thisChannel in range(self.__nChannels):
# bufferByChannel = self.__bufferList[thisChannel]
# #print self.__startIndex
# x = numpy.array(bufferByChannel)
# self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
#
# deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
# dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
# # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
#
# dataOut.data = self.__arrayBuffer
#
# self.__startIndex += self.__newNSamples
#
# return
class SSheightProfiles(Operation):
step = None
nsamples = None
bufferShape = None
profileShape = None
sshProfiles = None
profileIndex = None
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.isConfig = False
def setup(self,dataOut ,step = None , nsamples = None):
if step == None and nsamples == None:
raise ValueError("step or nheights should be specified ...")
self.step = step
self.nsamples = nsamples
self.__nChannels = dataOut.nChannels
self.__nProfiles = dataOut.nProfiles
self.__nHeis = dataOut.nHeights
shape = dataOut.data.shape #nchannels, nprofiles, nsamples
residue = (shape[1] - self.nsamples) % self.step
if residue != 0:
print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
numberProfile = self.nsamples
numberSamples = (shape[1] - self.nsamples)/self.step
self.bufferShape = int(shape[0]), int(numberSamples), int(numberProfile) # nchannels, nsamples , nprofiles
self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples) # nchannels, nprofiles, nsamples
self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
def run(self, dataOut, step, nsamples, code = None, repeat = None):
dataOut.flagNoData = True
profileIndex = None
#print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights)
#print(dataOut.getFreqRange(1)/1000.)
#exit(1)
if dataOut.flagDataAsBlock:
dataOut.data = numpy.average(dataOut.data,axis=1)
#print("jee")
dataOut.flagDataAsBlock = False
if not self.isConfig:
self.setup(dataOut, step=step , nsamples=nsamples)
#print("Setup done")
self.isConfig = True
if code is not None:
code = numpy.array(code)
code_block = code
if repeat is not None:
code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
#print(code_block.shape)
for i in range(self.buffer.shape[1]):
if code is not None:
self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block
else:
self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
#self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
for j in range(self.buffer.shape[0]):
self.sshProfiles[j] = numpy.transpose(self.buffer[j])
profileIndex = self.nsamples
deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
ippSeconds = (deltaHeight*1.0e-6)/(0.15)
#print("ippSeconds, dH: ",ippSeconds,deltaHeight)
try:
if dataOut.concat_m is not None:
ippSeconds= ippSeconds/float(dataOut.concat_m)
#print "Profile concat %d"%dataOut.concat_m
except:
pass
dataOut.data = self.sshProfiles
dataOut.flagNoData = False
dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
dataOut.profileIndex = profileIndex
dataOut.flagDataAsBlock = True
dataOut.ippSeconds = ippSeconds
dataOut.step = self.step
#print(numpy.shape(dataOut.data))
#exit(1)
#print("new data shape and time:", dataOut.data.shape, dataOut.utctime)
return dataOut
################################################################################3############################3
################################################################################3############################3
################################################################################3############################3
################################################################################3############################3
class SSheightProfiles2(Operation):
'''
Procesa por perfiles y por bloques
'''
bufferShape = None
profileShape = None
sshProfiles = None
profileIndex = None
#nsamples = None
#step = None
#deltaHeight = None
#init_range = None
__slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels',
'__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights')
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.isConfig = False
def setup(self,dataOut ,step = None , nsamples = None):
if step == None and nsamples == None:
raise ValueError("step or nheights should be specified ...")
self.step = step
self.nsamples = nsamples
self.__nChannels = int(dataOut.nChannels)
self.__nProfiles = int(dataOut.nProfiles)
self.__nHeis = int(dataOut.nHeights)
residue = (self.__nHeis - self.nsamples) % self.step
if residue != 0:
print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
self.init_range = dataOut.heightList[0]
#numberProfile = self.nsamples
numberSamples = (self.__nHeis - self.nsamples)/self.step
self.new_nHeights = numberSamples
self.bufferShape = int(self.__nChannels), int(numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples) # nchannels, nprofiles, nsamples
self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
def getNewProfiles(self, data, code=None, repeat=None):
if code is not None:
code = numpy.array(code)
code_block = code
if repeat is not None:
code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
if data.ndim == 2:
data = data.reshape(1,1,self.__nHeis )
#print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape)
for i in range(int(self.new_nHeights)): #nuevas alturas
if code is not None:
self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]*code_block
else:
self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
for j in range(self.__nChannels): #en los cananles
self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
#print("new profs Done")
def run(self, dataOut, step, nsamples, code = None, repeat = None):
if dataOut.flagNoData == True:
return dataOut
dataOut.flagNoData = True
#print("init data shape:", dataOut.data.shape)
#print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
# int(dataOut.nProfiles),int(dataOut.nHeights)))
profileIndex = None
# if not dataOut.flagDataAsBlock:
# dataOut.nProfiles = 1
if not self.isConfig:
self.setup(dataOut, step=step , nsamples=nsamples)
#print("Setup done")
self.isConfig = True
dataBlock = None
nprof = 1
if dataOut.flagDataAsBlock:
nprof = int(dataOut.nProfiles)
#print("dataOut nProfiles:", dataOut.nProfiles)
for profile in range(nprof):
if dataOut.flagDataAsBlock:
#print("read blocks")
self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
else:
#print("read profiles")
self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe
if profile == 0:
dataBlock = self.sshProfiles.copy()
else: #by blocks
dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
#print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
profileIndex = self.nsamples
#deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
ippSeconds = (self.deltaHeight*1.0e-6)/(0.15)
dataOut.data = dataBlock
#print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
dataOut.ippSeconds = ippSeconds
dataOut.step = self.step
dataOut.flagNoData = False
if dataOut.flagDataAsBlock:
dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
else:
dataOut.nProfiles = int(self.nsamples)
dataOut.profileIndex = dataOut.nProfiles
dataOut.flagDataAsBlock = True
dataBlock = None
#print("new data shape:", dataOut.data.shape, dataOut.utctime)
return dataOut
#import skimage.color
#import skimage.io
#import matplotlib.pyplot as plt
class removeProfileByFaradayHS(Operation):
'''
'''
__buffer_data = []
__buffer_times = []
buffer = None
outliers_IDs_list = []
__slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
'__dh','first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
'__count_exec','__initime','__dataReady','__ipp')
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.isConfig = False
def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3, minHei=None, maxHei=None):
if n == None and timeInterval == None:
raise ValueError("nprofiles or timeInterval should be specified ...")
if n != None:
self.n = n
self.navg = navg
self.profileMargin = profileMargin
self.thHistOutlier = thHistOutlier
self.__profIndex = 0
self.buffer = None
self._ipp = dataOut.ippSeconds
self.n_prof_released = 0
self.heightList = dataOut.heightList
self.init_prof = 0
self.end_prof = 0
self.__count_exec = 0
self.__profIndex = 0
self.first_utcBlock = None
self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
minHei = minHei
maxHei = maxHei
if minHei==None :
minHei = dataOut.heightList[0]
if maxHei==None :
maxHei = dataOut.heightList[-1]
self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
self.nChannels = dataOut.nChannels
self.nHeights = dataOut.nHeights
def filterSatsProfiles(self):
data = self.__buffer_data
#print(data.shape)
nChannels, profiles, heights = data.shape
indexes=[]
outliers_IDs=[]
for c in range(nChannels):
for h in range(self.minHei_idx, self.maxHei_idx):
power = data[c,:,h] * numpy.conjugate(data[c,:,h])
power = power.real
#power = (numpy.abs(data[c,:,h].real))
sortdata = numpy.sort(power, axis=None)
sortID=power.argsort()
index = _noise.hildebrand_sekhon2(sortdata,self.navg) #0.75-> buen valor
indexes.append(index)
outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
# print(outliers_IDs)
# fig,ax = plt.subplots()
# #ax.set_title(str(k)+" "+str(j))
# x=range(len(sortdata))
# ax.scatter(x,sortdata)
# ax.axvline(index)
# plt.grid()
# plt.show()
outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
outliers_IDs = numpy.unique(outliers_IDs)
outs_lines = numpy.sort(outliers_IDs)
# #print("outliers Ids: ", outs_lines, outs_lines.shape)
#hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
#Agrupando el histograma de outliers,
#my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
my_bins = numpy.linspace(0,9600, 96, endpoint=False)
hist, bins = numpy.histogram(outs_lines,bins=my_bins)
hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
#print(hist_outliers_indexes[0])
bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
#print(bins_outliers_indexes)
outlier_loc_index = []
# for n in range(len(bins_outliers_indexes)-1):
# for k in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin):
# outlier_loc_index.append(k)
outlier_loc_index = [e for n in range(len(bins_outliers_indexes)-1) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin) ]
outlier_loc_index = numpy.asarray(outlier_loc_index)
#print(len(numpy.unique(outlier_loc_index)), numpy.unique(outlier_loc_index))
# x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
# fig, ax = plt.subplots(1,2,figsize=(8, 6))
#
# dat = data[0,:,:].real
# m = numpy.nanmean(dat)
# o = numpy.nanstd(dat)
# #print(m, o, x.shape, y.shape)
# c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
# ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
# fig.colorbar(c)
# ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
# ax[1].hist(outs_lines,bins=my_bins)
# plt.show()
self.outliers_IDs_list = numpy.unique(outlier_loc_index)
return data
def cleanOutliersByBlock(self):
#print(self.__buffer_data[0].shape)
data = self.__buffer_data#.copy()
#print("cleaning shape inpt: ",data.shape)
'''
self.__buffer_data = []
spectrum = numpy.fft.fft2(data, axes=(0,2))
#print("spc : ",spectrum.shape)
(nch,nsamples, nh) = spectrum.shape
data2 = None
#print(data.shape)
spectrum2 = spectrum.copy()
for ch in range(nch):
dh = self.__dh
dt1 = (dh*1.0e-6)/(0.15)
dt2 = self.__buffer_times[1]-self.__buffer_times[0]
freqv = numpy.fft.fftfreq(nh, d=dt1)
freqh = numpy.fft.fftfreq(self.n, d=dt2)
#print("spc loop: ")
x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
z = numpy.abs(spectrum[ch,:,:])
# Find all peaks higher than the 98th percentile
peaks = z < numpy.percentile(z, 98)
#print(peaks)
# Set those peak coefficients to zero
spectrum2 = spectrum2 * peaks.astype(int)
data2 = numpy.fft.ifft2(spectrum2)
dat = numpy.log10(z.T)
dat2 = numpy.log10(spectrum2.T)
# m = numpy.mean(dat)
# o = numpy.std(dat)
# fig, ax = plt.subplots(2,1,figsize=(8, 6))
#
# c = ax[0].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
# #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
# date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
# #strftime('%Y-%m-%d %H:%M:%S')
# ax[0].set_title('Spectrum magnitude '+date_time)
# fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
#
#
# c = ax[1].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
# fig.colorbar(c)
# plt.show()
#print(data2.shape)
data = data2
#cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
'''
#print("cleanOutliersByBlock Done")
return self.filterSatsProfiles()
def fillBuffer(self, data, datatime):
if self.__profIndex == 0:
self.__buffer_data = data.copy()
else:
self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
self.__profIndex += 1
#self.__buffer_times.append(datatime)
def getData(self, data, datatime=None):
if self.__profIndex == 0:
self.__initime = datatime
self.__dataReady = False
self.fillBuffer(data, datatime)
dataBlock = None
if self.__profIndex == self.n:
#print("apnd : ",data)
#dataBlock = self.cleanOutliersByBlock()
dataBlock = self.filterSatsProfiles()
self.__dataReady = True
return dataBlock
if dataBlock is None:
return None, None
return dataBlock
def releaseBlock(self):
if self.n % self.lenProfileOut != 0:
raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
return None
data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
self.init_prof = self.end_prof
self.end_prof += self.lenProfileOut
#print("data release shape: ",dataOut.data.shape, self.end_prof)
self.n_prof_released += 1
#print("f_no_data ", dataOut.flagNoData)
return data
def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None):
#print("run op buffer 2D",dataOut.ippSeconds)
# self.nChannels = dataOut.nChannels
# self.nHeights = dataOut.nHeights
if not self.isConfig:
#print("init p idx: ", dataOut.profileIndex )
self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,
thHistOutlier=th_hist_outlier,minHei=minHei, maxHei=maxHei)
self.isConfig = True
dataBlock = None
if not dataOut.buffer_empty: #hay datos acumulados
if self.init_prof == 0:
self.n_prof_released = 0
self.lenProfileOut = nProfilesOut
dataOut.flagNoData = False
#print("tp 2 ",dataOut.data.shape)
self.init_prof = 0
self.end_prof = self.lenProfileOut
dataOut.nProfiles = self.lenProfileOut
if nProfilesOut == 1:
dataOut.flagDataAsBlock = False
else:
dataOut.flagDataAsBlock = True
#print("prof: ",self.init_prof)
dataOut.flagNoData = False
if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
#print("omitting: ", self.n_prof_released)
dataOut.flagNoData = True
dataOut.ippSeconds = self._ipp
dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
# print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
#dataOut.data = self.releaseBlock()
#########################################################3
if self.n % self.lenProfileOut != 0:
raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
return None
dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
self.init_prof = self.end_prof
self.end_prof += self.lenProfileOut
#print("data release shape: ",dataOut.data.shape, self.end_prof)
self.n_prof_released += 1
if self.end_prof >= (self.n +self.lenProfileOut):
self.init_prof = 0
self.__profIndex = 0
self.buffer = None
dataOut.buffer_empty = True
self.outliers_IDs_list = []
self.n_prof_released = 0
dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
#print("cleaning...", dataOut.buffer_empty)
dataOut.profileIndex = 0 #self.lenProfileOut
####################################################################
return dataOut
#print("tp 223 ",dataOut.data.shape)
dataOut.flagNoData = True
try:
#dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
self.__count_exec +=1
except Exception as e:
print("Error getting profiles data",self.__count_exec )
print(e)
sys.exit()
if self.__dataReady:
#print("omitting: ", len(self.outliers_IDs_list))
self.__count_exec = 0
#dataOut.data =
#self.buffer = numpy.flip(dataBlock, axis=1)
self.buffer = dataBlock
self.first_utcBlock = self.__initime
dataOut.utctime = self.__initime
dataOut.nProfiles = self.__profIndex
#dataOut.flagNoData = False
self.init_prof = 0
self.__profIndex = 0
self.__initime = None
dataBlock = None
self.__buffer_times = []
dataOut.error = False
dataOut.useInputBuffer = True
dataOut.buffer_empty = False
#print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
#print(self.__count_exec)
return dataOut
class RemoveProfileSats(Operation):
'''
Omite los perfiles contaminados con señal de satelites,
In: minHei = min_sat_range
max_sat_range
min_hei_ref
max_hei_ref
th = diference between profiles mean, ref and sats
Out:
profile clean
'''
isConfig = False
min_sats = 0
max_sats = 999999999
min_ref= 0
max_ref= 9999999999
needReshape = False
count = 0
thdB = 0
byRanges = False
min_sats = None
max_sats = None
noise = 0
def __init__(self, **kwargs):
Operation.__init__(self, **kwargs)
self.isConfig = False
def setup(self, dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList):
if rangeHeiList!=None:
self.byRanges = True
else:
if minHei==None or maxHei==None :
raise ValueError("Parameters heights are required")
if minRef==None or maxRef==None:
raise ValueError("Parameters heights are required")
if self.byRanges:
self.min_sats = []
self.max_sats = []
for min,max in rangeHeiList:
a,b = getHei_index(min, max, dataOut.heightList)
self.min_sats.append(a)
self.max_sats.append(b)
else:
self.min_sats, self.max_sats = getHei_index(minHei, maxHei, dataOut.heightList)
self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
self.th = th
self.thdB = thdB
self.isConfig = True
def compareRanges(self,data, minHei,maxHei):
# ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref])
# p_ref = 10*numpy.log10(ref.real)
# m_ref = numpy.mean(p_ref)
m_ref = self.noise
sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
p_sats = 10*numpy.log10(sats.real)
m_sats = numpy.mean(p_sats)
if m_sats > (m_ref + self.th): #and (m_sats > self.thdB):
#print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref))
#print("Removing profiles...")
return False
return True
def isProfileClean(self, data):
'''
Analiza solo 1 canal, y descarta todos...
'''
clean = True
if self.byRanges:
for n in range(len(self.min_sats)):
c = self.compareRanges(data,self.min_sats[n],self.max_sats[n])
clean = clean and c
else:
clean = (self.compareRanges(data, self.min_sats,self.max_sats))
return clean
def run(self, dataOut, minHei=None, maxHei=None, minRef=None, maxRef=None, th=5, thdB=65, rangeHeiList=None):
dataOut.flagNoData = True
if not self.isConfig:
self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
self.isConfig = True
#print(self.min_sats,self.max_sats)
if dataOut.flagDataAsBlock:
raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
else:
self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref))
if not self.isProfileClean(dataOut.data):
return dataOut
#dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)
#self.count += 1
dataOut.flagNoData = False
return dataOut