diff --git a/schainc/_noise.c b/schainc/_noise.c index 66ee196..38fedc9 100644 --- a/schainc/_noise.c +++ b/schainc/_noise.c @@ -57,7 +57,6 @@ static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) { if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) { return NULL; } - data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY); if (data_array == NULL) { @@ -89,18 +88,15 @@ static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) { j = j + 1; } - //double lnoise = sump / j; - Py_DECREF(data_array); return PyLong_FromLong(j); - } static PyMethodDef noiseMethods[] = { { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" }, - { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Get index for satellite cleaning" }, + { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Variation for satellite cleaning" }, { NULL, NULL, 0, NULL } }; diff --git a/schainpy/controller.py b/schainpy/controller.py index 087b7a2..9cbf8fe 100644 --- a/schainpy/controller.py +++ b/schainpy/controller.py @@ -266,10 +266,8 @@ class ProcUnitConf(ConfBase): def run(self): ''' ''' - return self.object.call(**self.getKwargs()) - class ReadUnitConf(ProcUnitConf): ELEMENTNAME = 'ReadUnit' @@ -633,18 +631,47 @@ class Project(Process): err = False n = len(self.configurations) - + flag_no_read = False + nProc_noRead = 0 while not err: #print("STAR") + n_proc = 0 for conf in self.getUnits(): #print("CONF: ",conf) - ok = conf.run() + #print( n_proc, nProc_noRead,flag_no_read) + if flag_no_read: + if n_proc >= nProc_noRead: + ok = conf.run() + else: + #print("break") + #print(ok, n_proc, nProc_noRead, n, flag_no_read) + n_proc += 1 + continue + else: + ok = conf.run() + + n_proc += 1 + #print(ok, n_proc, nProc_noRead, n, flag_no_read) + if ok == 'Error': #self.removeProcUnit(conf.id) #remove proc Unit n -= 1 continue + + elif ok == 'no_Read' and (not flag_no_read): + nProc_noRead = n_proc - 1 + flag_no_read = True + print("No read") + continue + elif ok == 'new_Read': + nProc_noRead = 0 + flag_no_read = False + print("read again") + continue elif not ok: + #print("not ok",ok) break + if n == 0: err = True diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 3f70377..fe08eae 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -159,6 +159,9 @@ class GenericData(object): class JROData(GenericData): + useInputBuffer = False + buffer_empty = True + systemHeaderObj = SystemHeader() radarControllerHeaderObj = RadarControllerHeader() type = None diff --git a/schainpy/model/graphics/jroplot_base.py b/schainpy/model/graphics/jroplot_base.py index 9c6f393..0f80196 100644 --- a/schainpy/model/graphics/jroplot_base.py +++ b/schainpy/model/graphics/jroplot_base.py @@ -434,6 +434,7 @@ class Plot(Operation): ax.firsttime = False if self.grid: ax.grid(True) + if not self.polar: ax.set_title('{} {} {}'.format( self.titles[n], @@ -442,6 +443,7 @@ class Plot(Operation): self.time_label), size=8) else: + ax.set_title('{}'.format(self.titles[n]), size=8) ax.set_ylim(0, 90) ax.set_yticks(numpy.arange(0, 90, 20)) @@ -480,6 +482,7 @@ class Plot(Operation): fig.canvas.manager.set_window_title('{} - {}'.format(self.title, self.getDateTime(self.data.max_time).strftime('%Y/%m/%d'))) + fig.canvas.draw() if self.show: fig.show() diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 6b0c4a7..2e260c3 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -11,7 +11,7 @@ import numpy from schainpy.model.graphics.jroplot_base import Plot, plt, log from itertools import combinations - +from matplotlib.ticker import LinearLocator class SpectraPlot(Plot): ''' @@ -86,7 +86,7 @@ class SpectraPlot(Plot): data = self.data[-1] z = data['spc'] - + print(z.shape, x.shape, y.shape) for n, ax in enumerate(self.axes): noise = data['noise'][n] if self.CODE == 'spc_moments': @@ -494,11 +494,13 @@ class SpectraCutPlot(Plot): self.height = 4.8 * self.nrows self.ylabel = 'Power [dB]' self.colorbar = False - self.plots_adjust.update({'left':0.15, 'hspace':0.3, 'right': 0.85, 'bottom':0.08}) + self.plots_adjust.update({'left':0.05, 'hspace':0.3, 'right': 0.9, 'bottom':0.08}) if len(self.selectedHeightsList) > 0: self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight)) + + def update(self, dataOut): if len(self.channelList) == 0: self.channelList = dataOut.channelList @@ -544,21 +546,33 @@ class SpectraCutPlot(Plot): else: index = numpy.arange(0, len(y), int((len(y))/9)) #print("inde x ", index, self.axes) + for n, ax in enumerate(self.axes): + if ax.firsttime: + + self.xmax = self.xmax if self.xmax else numpy.nanmax(x) self.xmin = self.xmin if self.xmin else -self.xmax self.ymin = self.ymin if self.ymin else numpy.nanmin(z) self.ymax = self.ymax if self.ymax else numpy.nanmax(z) + + ax.plt = ax.plot(x, z[n, :, index].T) labels = ['Range = {:2.1f}km'.format(y[i]) for i in index] self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8}) + ax.minorticks_on() + ax.grid(which='major', axis='both') + ax.grid(which='minor', axis='x') else: for i, line in enumerate(ax.plt): line.set_data(x, z[n, :, index[i]]) + + self.titles.append('CH {}'.format(self.channelList[n])) plt.suptitle(self.maintitle, fontsize=10) + class BeaconPhase(Plot): __isConfig = None @@ -937,7 +951,7 @@ class NoiselessRTIPlot(Plot): if nch != 1: aux = [] for c in self.channelList: - aux.append(n0[c]) + aux.append(n0[c]) n0 = numpy.asarray(aux) noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh)) #print(dataOut.elevationList, dataOut.azimuthList) diff --git a/schainpy/model/io/utils.py b/schainpy/model/io/utils.py index 2660a31..d1ea2a0 100644 --- a/schainpy/model/io/utils.py +++ b/schainpy/model/io/utils.py @@ -4,6 +4,7 @@ Utilities for IO modules import os from datetime import datetime +import numpy def folder_in_range(folder, start_date, end_date, pattern): """ @@ -22,3 +23,29 @@ def folder_in_range(folder, start_date, end_date, pattern): except: raise ValueError('Folder {} does not match {} format'.format(folder, pattern)) return start_date <= dt.date() <= end_date + + +def getHei_index( minHei, maxHei, heightList): + if (minHei < heightList[0]): + minHei = heightList[0] + + if (maxHei > heightList[-1]): + maxHei = heightList[-1] + + minIndex = 0 + maxIndex = 0 + heights = numpy.asarray(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(heightList) + return minIndex,maxIndex diff --git a/schainpy/model/proc/jroproc_base.py b/schainpy/model/proc/jroproc_base.py index 3a7ebba..4390e56 100644 --- a/schainpy/model/proc/jroproc_base.py +++ b/schainpy/model/proc/jroproc_base.py @@ -53,18 +53,24 @@ class ProcessingUnit(object): return self.operations[objId] def call(self, **kwargs): - ''' - ''' + mybool = (self.dataOut.type == 'Voltage') and self.dataOut.useInputBuffer and (not self.dataOut.buffer_empty) #liberar desde buffer try: - if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error: - return self.dataIn.isReady() - elif self.dataIn is None or not self.dataIn.error: + + if mybool: + #print("run jeje") self.run(**kwargs) - elif self.dataIn.error: - self.dataOut.error = self.dataIn.error - self.dataOut.flagNoData = True + else: + if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error: + return self.dataIn.isReady() + elif self.dataIn is None or not self.dataIn.error: #unidad de lectura o procesamiento regular + self.run(**kwargs) + elif self.dataIn.error: + self.dataOut.error = self.dataIn.error + self.dataOut.flagNoData = True + print("exec proc error") + except: err = traceback.format_exc() @@ -76,8 +82,10 @@ class ProcessingUnit(object): log.error(err, self.name) self.dataOut.error = True + for op, optype, opkwargs in self.operations: - if optype == 'other' and self.dataOut.isReady(): + + if (optype == 'other' and self.dataOut.isReady()) or mybool: try: self.dataOut = op.run(self.dataOut, **opkwargs) except Exception as e: @@ -89,7 +97,21 @@ class ProcessingUnit(object): elif optype == 'external' and self.dataOut.error: op.queue.put(copy.deepcopy(self.dataOut)) - return 'Error' if self.dataOut.error else self.dataOut.isReady() + + if not self.dataOut.error: + if self.dataOut.type == 'Voltage': + if not self.dataOut.buffer_empty : #continue + return 'no_Read' + elif self.dataOut.useInputBuffer and (self.dataOut.buffer_empty) and self.dataOut.isReady() : + return 'new_Read' + else: + return True + else: + #print("ret True") + return True + else: + return 'Error' + #return 'Error' if self.dataOut.error else True #self.dataOut.isReady() def setup(self): diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 2a557ac..380c4c6 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -39,6 +39,8 @@ class SpectraProc(ProcessingUnit): def __updateSpecFromVoltage(self): + + self.dataOut.timeZone = self.dataIn.timeZone self.dataOut.dstFlag = self.dataIn.dstFlag self.dataOut.errorCount = self.dataIn.errorCount @@ -130,7 +132,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.copy(self.dataIn) except Exception as e: - print(e) + print("Error dataIn ",e) if shift_fft: #desplaza a la derecha en el eje 2 determinadas posiciones @@ -158,7 +160,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.ippFactor = 1 self.dataOut.nFFTPoints = nFFTPoints - + #print(" volts ch,prof, h: ", self.dataIn.data.shape) if self.buffer is None: self.buffer = numpy.zeros((self.dataIn.nChannels, nProfiles, @@ -195,6 +197,7 @@ class SpectraProc(ProcessingUnit): self.firstdatatime = self.dataIn.utctime if self.profIndex == nProfiles: + self.__updateSpecFromVoltage() if pairsList == None: self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)] @@ -209,6 +212,7 @@ class SpectraProc(ProcessingUnit): raise ValueError("The type of input object '%s' is not valid".format( self.dataIn.type)) + def __selectPairs(self, pairsList): if not pairsList: @@ -492,14 +496,14 @@ class removeDC(Operation): return self.dataOut class getNoise(Operation): - + warnings = False def __init__(self): Operation.__init__(self) - def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None): + def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False): self.dataOut = dataOut - + self.warnings = warnings if minHei == None: minHei = self.dataOut.heightList[0] @@ -507,13 +511,15 @@ class getNoise(Operation): maxHei = self.dataOut.heightList[-1] if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei): - print('minHei: %.2f is out of the heights range' % (minHei)) - print('minHei is setting to %.2f' % (self.dataOut.heightList[0])) + if self.warnings: + print('minHei: %.2f is out of the heights range' % (minHei)) + print('minHei is setting to %.2f' % (self.dataOut.heightList[0])) minHei = self.dataOut.heightList[0] if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei): - print('maxHei: %.2f is out of the heights range' % (maxHei)) - print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1])) + if self.warnings: + print('maxHei: %.2f is out of the heights range' % (maxHei)) + print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1])) maxHei = self.dataOut.heightList[-1] @@ -535,13 +541,15 @@ class getNoise(Operation): maxFreq = freqrange[-1] if (minFreq < freqrange[0]) or (minFreq > maxFreq): - print('minFreq: %.2f is out of the frequency range' % (minFreq)) - print('minFreq is setting to %.2f' % (freqrange[0])) + if self.warnings: + print('minFreq: %.2f is out of the frequency range' % (minFreq)) + print('minFreq is setting to %.2f' % (freqrange[0])) minFreq = freqrange[0] if (maxFreq > freqrange[-1]) or (maxFreq < minFreq): - print('maxFreq: %.2f is out of the frequency range' % (maxFreq)) - print('maxFreq is setting to %.2f' % (freqrange[-1])) + if self.warnings: + print('maxFreq: %.2f is out of the frequency range' % (maxFreq)) + print('maxFreq is setting to %.2f' % (freqrange[-1])) maxFreq = freqrange[-1] indminPoint = numpy.where(freqrange >= minFreq) @@ -557,13 +565,15 @@ class getNoise(Operation): maxVel = velrange[-1] if (minVel < velrange[0]) or (minVel > maxVel): - print('minVel: %.2f is out of the velocity range' % (minVel)) - print('minVel is setting to %.2f' % (velrange[0])) + if self.warnings: + print('minVel: %.2f is out of the velocity range' % (minVel)) + print('minVel is setting to %.2f' % (velrange[0])) minVel = velrange[0] if (maxVel > velrange[-1]) or (maxVel < minVel): - print('maxVel: %.2f is out of the velocity range' % (maxVel)) - print('maxVel is setting to %.2f' % (velrange[-1])) + if self.warnings: + print('maxVel: %.2f is out of the velocity range' % (maxVel)) + print('maxVel is setting to %.2f' % (velrange[-1])) maxVel = velrange[-1] indminPoint = numpy.where(velrange >= minVel) @@ -1094,14 +1104,14 @@ class IntegrationFaradaySpectra(Operation): n = None minHei_ind = None maxHei_ind = None - avg = 1.0 + navg = 1.0 factor = 0.0 def __init__(self): Operation.__init__(self) - def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, factor=0.75): + def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75): """ Set the parameters of the integration class. @@ -1125,6 +1135,7 @@ class IntegrationFaradaySpectra(Operation): self.__byTime = False self.factor = factor + self.navg = avg #self.ByLags = dataOut.ByLags ###REDEFINIR self.ByLags = False @@ -1277,13 +1288,11 @@ class IntegrationFaradaySpectra(Operation): outliers_IDs=outliers_IDs.astype(numpy.dtype('int64')) indexes=numpy.array(indexes) indexmin=numpy.min(indexes) - #print("clean CH: ", i) + if indexmin != buffer1.shape[0]: if self.nChannels > 1: cspc_outliers_exist= True - #print("outliers cspc") - ###sortdata=numpy.sort(buffer1,axis=0) - ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0) + lt=outliers_IDs avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0) @@ -1400,7 +1409,7 @@ class IntegrationFaradaySpectra(Operation): return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc - def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, factor=0.75): + def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75): self.dataOut = dataOut.copy() if n == 1: return self.dataOut @@ -1410,7 +1419,7 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc) if not self.isConfig: - self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, factor) + self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor) self.isConfig = True if not self.ByLags: @@ -1742,7 +1751,9 @@ class IncohInt(Operation): Add a profile to the __buffer_spc and increase in one the __profileIndex """ - + if data_spc.all() == numpy.nan : + print("nan ") + return self.__buffer_spc += data_spc if data_cspc is None: diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 7cd3b7a..0e821a4 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -4,9 +4,11 @@ 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 class VoltageProc(ProcessingUnit): @@ -19,12 +21,15 @@ class VoltageProc(ProcessingUnit): self.setupReq = False def run(self): - + #print("running volt proc") if self.dataIn.type == 'AMISR': self.__updateObjFromAmisrInput() - if self.dataIn.type == 'Voltage': - self.dataOut.copy(self.dataIn) + if self.dataOut.buffer_empty: + if self.dataIn.type == 'Voltage': + self.dataOut.copy(self.dataIn) + #print("new volts reading") + def __updateObjFromAmisrInput(self): @@ -1680,6 +1685,7 @@ class SSheightProfiles(Operation): dataOut.flagNoData = True profileIndex = None + #print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights) #print(dataOut.getFreqRange(1)/1000.) #exit(1) if dataOut.flagDataAsBlock: @@ -1691,30 +1697,11 @@ class SSheightProfiles(Operation): #print("Setup done") self.isConfig = True - #DC_Hae = numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j]) - #DC_Hae = numpy.array([ 0.001025 +0.0516375j, 0.03485 +0.20923125j, -0.168 -0.02720625j, - #-0.1105375 +0.0707125j, -0.20309375-0.09670625j, 0.189775 +0.02716875j])*(-3.5) - - #DC_Hae = numpy.array([ -32.26 +8.66j, -32.26 +8.66j]) - #DC_Hae = numpy.array([-2.78500000e-01 -1.39175j, -6.63237294e+02+210.4268625j]) - - #print(dataOut.data[0,13:15]) - #dataOut.data = dataOut.data - DC_Hae[:,None] - #print(dataOut.data[0,13:15]) - #exit(1) if code is not None: code = numpy.array(code) code_block = code - ''' - roll = 0 - code = numpy.roll(code,roll,axis=0) - code = numpy.reshape(code,(5,100,64)) - block = dataOut.CurrentBlock%5 - - day_dif = 0 #day_19_Oct_2021: 3 - code_block = code[block-1-day_dif,:,:] - ''' + if repeat is not None: code_block = numpy.repeat(code_block, repeats=repeat, axis=1) #print(code_block.shape) @@ -1756,3 +1743,492 @@ class SSheightProfiles(Operation): #exit(1) return dataOut +################################################################################3############################3 +################################################################################3############################3 +################################################################################3############################3 +################################################################################3############################3 + +class SSheightProfiles2(Operation): + ''' + Procesa por perfiles y por bloques + ''' + + 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 = 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)) + + deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + #numberProfile = self.nsamples + self.numberSamples = (self.__nHeis - self.nsamples)/self.step + + self.new_nHeights = self.numberSamples + + self.bufferShape = int(self.__nChannels), int(self.numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles + self.profileShape = int(self.__nChannels), int(self.nsamples), int(self.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) + #print("buff, data, :",self.buffer.shape, data.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]) + + + + + def run(self, dataOut, step, nsamples, code = None, repeat = None): + + 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: + self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat) + else: + self.getNewProfiles(dataOut.data, code=code, repeat=repeat) + 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 = (deltaHeight*1.0e-6)/(0.15) + + + dataOut.data = dataBlock + dataOut.heightList = numpy.arange(self.new_nHeights) *self.step*deltaHeight + dataOut.heightList[0] + #print("show me: ",dataOut.nHeights, self.new_nHeights) + #dataOut.nHeights = int(self.new_nHeights) + 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) + + return dataOut + + + + +#import skimage.color +#import skimage.io +import matplotlib.pyplot as plt + +class clean2DArray(Operation): + ''' + + ''' + isConfig = False + n = None + __timeInterval = None + __profIndex = 0 + __byTime = False + __dataReady = False + __buffer_data = [] + __buffer_times = [] + __initime = None + __count_exec = 0 + __profIndex = 0 + buffer = None + lenProfileOut = 1 + + init_prof = 0 + end_prof = 0 + n_profiles = 0 + first_utcBlock = None + __dh = 0 + + def __init__(self, **kwargs): + + Operation.__init__(self, **kwargs) + self.isConfig = False + + def setup(self,dataOut, n=None , timeInterval=None): + + if n == None and timeInterval == None: + raise ValueError("nprofiles or timeInterval should be specified ...") + + if n != None: + self.n = n + self.__byTime = False + else: + self.__timeInterval = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line + self.n = 9999 + self.__byTime = True + + self.__profIndex = 0 + self.buffer = None + self.lenProfileOut = 1 + + self.init_prof = 0 + self.end_prof = 0 + self.n_profiles = 0 + self.first_utcBlock = None + self.__dh = dataOut.heightList[1] - dataOut.heightList[0] + + 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 + + #nch = + + 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,:,:]) + + + dat = numpy.log10(z.T) + m = numpy.mean(dat) + o = numpy.std(dat) + fig, ax = plt.subplots(figsize=(8, 6)) + + c = ax.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.set_title('Spectrum magnitude '+date_time) + fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time) + fig.colorbar(c) + + plt.show() + + + + #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape() + + print("cleanOutliersByBlock Done") + return data + + def byTime(self, data, datatime): + + self.__dataReady = False + + self.fillBuffer(data, datatime) + dataBlock = None + + if (datatime - self.__initime) >= self.__timeInterval: + dataBlock = self.cleanOutliersByBlock() + self.__dataReady = True + self.n = self.__profIndex + return dataBlock + + def byProfiles(self, data, datatime): + self.__dataReady = False + + self.fillBuffer(data, datatime) + dataBlock = None + + if self.__profIndex == self.n: + #print("apnd : ",data) + dataBlock = self.cleanOutliersByBlock() + self.__dataReady = True + + return dataBlock + + 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 + + if self.__byTime: + dataBlock = self.byTime(data, datatime) + else: + dataBlock = self.byProfiles(data, datatime) + + + + if dataBlock is None: + return None + + + + return dataBlock + + def releaseBlock(self, dataOut): + + if self.n % self.lenProfileOut != 0: + raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles)) + return None + + dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt + dataOut.profileIndex = self.lenProfileOut + dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds + self.init_prof = self.end_prof + self.end_prof += self.lenProfileOut + print("data release shape: ",dataOut.data.shape, self.end_prof) + + if self.end_prof >= (self.n +self.lenProfileOut): + self.init_prof = 0 + self.__profIndex = 0 + self.buffer = None + dataOut.buffer_empty = True + return dataOut + + def run(self, dataOut, n=None, timeInterval=None, nProfilesOut=1): + #print("run op buffer 2D") + self.nChannels = dataOut.nChannels + self.nHeights = dataOut.nHeights + + if not self.isConfig: + self.setup(dataOut,n=n, timeInterval=timeInterval) + self.isConfig = True + + dataBlock = None + + if not dataOut.buffer_empty: #hay datos acumulados + + if self.init_prof == 0: + self.lenProfileOut = nProfilesOut + dataOut.flagNoData = False + #print("tp 2 ",dataOut.data.shape) + #(ch, self.n_profiles, nh) = self.buffer.shape + #print("tp 3 ",self.dataOut.data.shape) + #print("rel: ",self.buffer[:,-1,:]) + self.init_prof = 0 + self.end_prof = self.lenProfileOut + self.first_utcBlock = dataOut.utctime + dataOut.nProfiles = self.lenProfileOut + dataOut.error = False + dataOut.flagDataAsBlock = True + #print("prof: ",self.init_prof) + dataOut.flagNoData = False + return self.releaseBlock(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: + self.__count_exec = 0 + #dataOut.data = + self.buffer = numpy.flip(dataBlock, axis=1) + + dataOut.utctime = self.__initime + dataOut.nProfiles = self.__profIndex + #dataOut.flagNoData = False + 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))) + #return None + + + #print(self.__count_exec) + + return dataOut + +class CleanProfileSats(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 + + 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) + + 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," mRef: ", m_ref, (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 + + if dataOut.flagDataAsBlock: + raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False") + + else: + + 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