diff --git a/.gitignore b/.gitignore index c78f88d..e0bc00e 100644 --- a/.gitignore +++ b/.gitignore @@ -112,5 +112,3 @@ schainpy/scripts/ .vscode trash *.log -schainpy/scripts/testDigitalRF.py -schainpy/scripts/testDigitalRFWriter.py diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 39844e7..bebece8 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -190,8 +190,7 @@ class JROData(GenericData): profileIndex = None error = None data = None - data_plt = None - + nmodes = None def __str__(self): @@ -462,6 +461,7 @@ class Spectra(JROData): ippFactor = None profileIndex = 0 plotting = "spectra" + def __init__(self): ''' Constructor @@ -554,9 +554,12 @@ class Spectra(JROData): deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor) velrange = deltav * (numpy.arange(self.nFFTPoints + - extrapoints) - self.nFFTPoints / 2.) # - deltav/2 + extrapoints) - self.nFFTPoints / 2.) - return velrange + if self.nmodes: + return velrange/self.nmodes + else: + return velrange def getNPairs(self): diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index d049627..b6d0046 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -7,14 +7,190 @@ from .plotting_codes import * from schainpy.model.proc.jroproc_base import MPDecorator from schainpy.utils import log -class FitGauPlot_(Figure): +class ParamLine_(Figure): + + isConfig = None + + def __init__(self): + + self.isConfig = False + self.WIDTH = 300 + self.HEIGHT = 200 + self.counter_imagwr = 0 + + def getSubplots(self): + + nrow = self.nplots + ncol = 3 + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + self.nplots = nplots + + self.createFigure(id=id, + wintitle=wintitle, + show=show) + + nrow,ncol = self.getSubplots() + colspan = 3 + rowspan = 1 + + for i in range(nplots): + self.addAxes(nrow, ncol, i, 0, colspan, rowspan) + + def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + yreal = y[channelIndexList,:].real + yimag = y[channelIndexList,:].imag + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity - IQ" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag)) + if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag)) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + + axes.pline(x, yreal[i,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2) + + def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) + yreal = y.real + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(yreal) + if ymax == None: ymax = numpy.nanmax(yreal) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = yreal[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, save=False, + figpath='./', figfile=None, show=True, wr_period=1, + ftp=False, server=None, folder=None, username=None, password=None): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + """ + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError("Channel %d is not in dataOut.channelList" % channel) + channelIndexList.append(dataOut.channelList.index(channel)) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) + + y = dataOut.RR + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = y[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + self.draw() + + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex) + figfile = self.getFilename(name = str_datetime) + + self.save(figpath=figpath, + figfile=figfile, + save=save, + ftp=ftp, + wr_period=wr_period, + thisDatetime=thisDatetime) + + + +class SpcParamPlot_(Figure): isConfig = None __nsubplots = None WIDTHPROF = None HEIGHTPROF = None - PREFIX = 'fitgau' + PREFIX = 'SpcParam' def __init__(self, **kwargs): Figure.__init__(self, **kwargs) @@ -83,7 +259,7 @@ class FitGauPlot_(Figure): save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False, - xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1): + xaxis="frequency", colormap='jet', normFactor=None , Selector = 0): """ @@ -119,23 +295,22 @@ class FitGauPlot_(Figure): # else: # factor = normFactor if xaxis == "frequency": - x = dataOut.spc_range[0] + x = dataOut.spcparam_range[0] xlabel = "Frequency (kHz)" elif xaxis == "time": - x = dataOut.spc_range[1] + x = dataOut.spcparam_range[1] xlabel = "Time (ms)" - else: - x = dataOut.spc_range[2] + else: + x = dataOut.spcparam_range[2] xlabel = "Velocity (m/s)" - ylabel = "Range (Km)" + ylabel = "Range (km)" y = dataOut.getHeiRange() - z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor - print('GausSPC', z[0,32,10:40]) + z = dataOut.SPCparam[Selector] /1966080.0#/ dataOut.normFactor#GauSelector] #dataOut.data_spc/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) zdB = 10*numpy.log10(z) @@ -657,7 +832,7 @@ class WindProfilerPlot_(Figure): # tmax = None x = dataOut.getTimeRange1(dataOut.paramInterval) - y = dataOut.heightList + y = dataOut.heightList z = dataOut.data_output.copy() nplots = z.shape[0] #Number of wind dimensions estimated nplotsw = nplots @@ -666,13 +841,14 @@ class WindProfilerPlot_(Figure): #If there is a SNR function defined if dataOut.data_SNR is not None: nplots += 1 - SNR = dataOut.data_SNR - SNRavg = numpy.average(SNR, axis=0) + SNR = dataOut.data_SNR[0] + SNRavg = SNR#numpy.average(SNR, axis=0) SNRdB = 10*numpy.log10(SNR) SNRavgdB = 10*numpy.log10(SNRavg) - if SNRthresh == None: SNRthresh = -5.0 + if SNRthresh == None: + SNRthresh = -5.0 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0] for i in range(nplotsw): @@ -741,8 +917,7 @@ class WindProfilerPlot_(Figure): axes = self.axesList[i*self.__nsubplots] z1 = z[i,:].reshape((1,-1))*windFactor[i] - #z1=numpy.ma.masked_where(z1==0.,z1) - + axes.pcolorbuffer(x, y, z1, xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i], xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, @@ -792,8 +967,8 @@ class ParametersPlot_(Figure): self.isConfig = False self.__nsubplots = 1 - self.WIDTH = 800 - self.HEIGHT = 180 + self.WIDTH = 300 + self.HEIGHT = 550 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 @@ -905,7 +1080,7 @@ class ParametersPlot_(Figure): # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) xlabel = "" - ylabel = "Range (Km)" + ylabel = "Range (km)" update_figfile = False @@ -949,24 +1124,81 @@ class ParametersPlot_(Figure): self.setWinTitle(title) - for i in range(self.nchan): - index = channelIndexList[i] - title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact] - z1 = z[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, z1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap=colormap) - - if showSNR: - title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact + 1] - SNRdB1 = SNRdB[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, SNRdB1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap='jet') +# for i in range(self.nchan): +# index = channelIndexList[i] +# title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact] +# z1 = z[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, z1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap=colormap) +# +# if showSNR: +# title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact + 1] +# SNRdB1 = SNRdB[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, SNRdB1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=0 + index = channelIndexList[i] + title = "Factor de reflectividad Z [dBZ]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap=colormap) + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=1 + index = channelIndexList[i] + title = "Velocidad vertical Doppler [m/s]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=2 + index = channelIndexList[i] + title = "Intensidad de lluvia [mm/h]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') self.draw() @@ -1067,9 +1299,8 @@ class Parameters1Plot_(Figure): save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, server=None, folder=None, username=None, password=None, ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): - #print inspect.getargspec(self.run).args - """ + """ Input: dataOut : id : diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 47fd698..1f6e1f3 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -42,6 +42,8 @@ class SpectraPlot_(Figure): self.__xfilter_ena = False self.__yfilter_ena = False + + self.indice=1 def getSubplots(self): @@ -139,10 +141,9 @@ class SpectraPlot_(Figure): x = dataOut.getVelRange(1) xlabel = "Velocity (m/s)" - ylabel = "Range (Km)" + ylabel = "Range (km)" y = dataOut.getHeiRange() - z = dataOut.data_spc/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) zdB = 10*numpy.log10(z) @@ -155,6 +156,7 @@ class SpectraPlot_(Figure): thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Spectra" + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith) @@ -223,6 +225,7 @@ class SpectraPlot_(Figure): ftp=ftp, wr_period=wr_period, thisDatetime=thisDatetime) + return dataOut @MPDecorator @@ -252,6 +255,8 @@ class CrossSpectraPlot_(Figure): self.EXP_CODE = None self.SUB_EXP_CODE = None self.PLOT_POS = None + + self.indice=0 def getSubplots(self): @@ -396,6 +401,7 @@ class CrossSpectraPlot_(Figure): self.isConfig = True self.setWinTitle(title) + for i in range(self.nplots): pair = dataOut.pairsList[pairsIndexList[i]] @@ -420,7 +426,7 @@ class CrossSpectraPlot_(Figure): xlabel=xlabel, ylabel=ylabel, title=title, ticksize=9, colormap=power_cmap, cblabel='') - coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:]) + coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] ) coherence = numpy.abs(coherenceComplex) # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi @@ -439,8 +445,6 @@ class CrossSpectraPlot_(Figure): xlabel=xlabel, ylabel=ylabel, title=title, ticksize=9, colormap=phase_cmap, cblabel='') - - self.draw() self.save(figpath=figpath, @@ -470,7 +474,7 @@ class RTIPlot_(Figure): self.__nsubplots = 1 self.WIDTH = 800 - self.HEIGHT = 180 + self.HEIGHT = 250 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 @@ -1497,9 +1501,6 @@ class BeaconPhase_(Figure): avgcoherenceComplex = ccf/numpy.sqrt(powa*powb) phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi - #print "Phase %d%d" %(pair[0], pair[1]) - #print phase[dataOut.beacon_heiIndexList] - if dataOut.beacon_heiIndexList: phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList]) else: diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 536daa2..9b4b1b6 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -434,7 +434,6 @@ def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='' def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''): ax = iplot.axes - printLabels(ax, xlabel, ylabel, title) for i in range(len(ax.lines)): diff --git a/schainpy/model/io/bltrIO_param.py b/schainpy/model/io/bltrIO_param.py index b6b1dc3..70a996c 100644 --- a/schainpy/model/io/bltrIO_param.py +++ b/schainpy/model/io/bltrIO_param.py @@ -111,7 +111,6 @@ class BLTRParamReader(JRODataReader, ProcessingUnit): timezone=0, status_value=0, **kwargs): - self.path = path self.startDate = startDate self.endDate = endDate diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index 2a6c9f4..941ea92 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -1815,7 +1815,7 @@ class JRODataWriter(JRODataIO): return 1 - def run(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs): + def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs): if not(self.isConfig): diff --git a/schainpy/model/io/jroIO_bltr.py b/schainpy/model/io/jroIO_bltr.py index dd29481..593d3cb 100644 --- a/schainpy/model/io/jroIO_bltr.py +++ b/schainpy/model/io/jroIO_bltr.py @@ -63,9 +63,6 @@ class Header(object): if attr: message += "%s = %s" % ("size", attr) + "\n" - # print message - - FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes ('FileMgcNumber', ' endFp: sys.stderr.write( "Warning %s: Size value read from System Header is lower than it has to be\n" % fp) @@ -537,9 +440,6 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa FileList.append(IndexFile) nFiles += 1 - # print 'Files2Read' - # print 'Existen '+str(nFiles)+' archivos .fdt' - self.filenameList = FileList # List of files from least to largest by names def run(self, **kwargs): @@ -553,7 +453,6 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.isConfig = True self.getData() - # print 'running' def setup(self, path=None, startDate=None, @@ -590,22 +489,19 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa if self.flagNoMoreFiles: self.dataOut.flagNoData = True - print('NoData se vuelve true') return 0 self.fp = self.path self.Files2Read(self.fp) self.readFile(self.fp) self.dataOut.data_spc = self.data_spc - self.dataOut.data_cspc = self.data_cspc - self.dataOut.data_output = self.data_output - - print('self.dataOut.data_output', shape(self.dataOut.data_output)) - - # self.removeDC() - return self.dataOut.data_spc - - def readFile(self, fp): + self.dataOut.data_cspc =self.data_cspc + self.dataOut.data_output=self.data_output + + return self.dataOut.data_spc + + + def readFile(self,fp): ''' You must indicate if you are reading in Online or Offline mode and load the The parameters for this file reading mode. @@ -615,23 +511,18 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa 1. Get the BLTR FileHeader. 2. Start reading the first block. ''' - - # The address of the folder is generated the name of the .fdt file that will be read - print("File: ", self.fileSelector + 1) - + if self.fileSelector < len(self.filenameList): self.fpFile = str(fp) + '/' + \ str(self.filenameList[self.fileSelector]) - # print self.fpFile fheader = FileHeaderBLTR() fheader.FHread(self.fpFile) # Bltr FileHeader Reading self.nFDTdataRecors = fheader.nFDTdataRecors self.readBlock() # Block reading else: - print('readFile FlagNoData becomes true') - self.flagNoMoreFiles = True + self.flagNoMoreFiles=True self.dataOut.flagNoData = True return 0 @@ -658,12 +549,11 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa 2. Fill the buffer with the current block number. ''' - - if self.BlockCounter < self.nFDTdataRecors - 2: - print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') - if self.ReadMode == 1: - rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1) - elif self.ReadMode == 0: + + if self.BlockCounter < self.nFDTdataRecors-1: + if self.ReadMode==1: + rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1) + elif self.ReadMode==0: rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter) rheader.RHread(self.fpFile) # Bltr FileHeader Reading @@ -683,31 +573,26 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.nRdPairs = len(self.dataOut.pairsList) self.dataOut.nRdPairs = self.nRdPairs - - self.__firstHeigth = rheader.StartRangeSamp - self.__deltaHeigth = rheader.SampResolution - self.dataOut.heightList = self.__firstHeigth + \ - numpy.array(list(range(self.nHeights))) * self.__deltaHeigth - self.dataOut.channelList = list(range(self.nChannels)) - self.dataOut.nProfiles = rheader.nProfiles - self.dataOut.nIncohInt = rheader.nIncohInt - self.dataOut.nCohInt = rheader.nCohInt - self.dataOut.ippSeconds = 1 / float(rheader.PRFhz) - self.dataOut.PRF = rheader.PRFhz - self.dataOut.nFFTPoints = rheader.nProfiles - self.dataOut.utctime = rheader.nUtime - self.dataOut.timeZone = 0 - self.dataOut.normFactor = self.dataOut.nProfiles * \ - self.dataOut.nIncohInt * self.dataOut.nCohInt - self.dataOut.outputInterval = self.dataOut.ippSeconds * \ - self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles - - self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN - print('self.data_output', shape(self.data_output)) - self.dataOut.velocityX = [] - self.dataOut.velocityY = [] - self.dataOut.velocityV = [] - + self.__firstHeigth=rheader.StartRangeSamp + self.__deltaHeigth=rheader.SampResolution + self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth + self.dataOut.channelList = range(self.nChannels) + self.dataOut.nProfiles=rheader.nProfiles + self.dataOut.nIncohInt=rheader.nIncohInt + self.dataOut.nCohInt=rheader.nCohInt + self.dataOut.ippSeconds= 1/float(rheader.PRFhz) + self.dataOut.PRF=rheader.PRFhz + self.dataOut.nFFTPoints=rheader.nProfiles + self.dataOut.utctime=rheader.nUtime + self.dataOut.timeZone=0 + self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt + self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles + + self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN + self.dataOut.velocityX=[] + self.dataOut.velocityY=[] + self.dataOut.velocityV=[] + '''Block Reading, the Block Data is received and Reshape is used to give it shape. ''' @@ -734,18 +619,17 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa y = rho * numpy.sin(phi) return(x, y) - if self.DualModeIndex == self.ReadMode: - - self.data_fft = numpy.fromfile( - startDATA, [('complex', ' 0.0001) : -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) -# -# if numpy.amax(popt)>numpy.amax(yMean)*0.3: -# FitGauss=gaus(xSamples,*popt) -# -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Dentro', Height -# except RuntimeError: -# -# try: -# for j in range(len(ySamples[1])): -# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]])) -# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion1', Height -# except RuntimeError: -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion2', Height -# except RuntimeError: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Exepcion3', Height -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# #print 'Verificador: Fuera', Height -# -# -# -# Maximun=numpy.amax(yMean) -# eMinus1=Maximun*numpy.exp(-1) -# -# HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) -# HalfWidth= xFrec[HWpos] -# GCpos=Find(FitGauss, numpy.amax(FitGauss)) -# Vpos=Find(FactNorm, numpy.amax(FactNorm)) -# #Vpos=numpy.sum(FactNorm)/len(FactNorm) -# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) ))) -# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos -# '''****** Getting Fij ******''' -# -# GaussCenter=xFrec[GCpos] -# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): -# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 -# else: -# Fij=abs(GaussCenter-HalfWidth)+0.0000001 -# -# '''****** Getting Frecuency range of significant data ******''' -# -# Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) -# -# if Rangpos5 and len(FrecRange) 0.: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag -# #print 'Vmag',Vmag -# else: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN) -# -# if abs(Vx)<100 and abs(Vx) > 0.: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang -# #print 'Vang',Vang -# else: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN) -# -# if abs(GaussCenter)<2: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos]) -# -# else: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN) -# -# -# # print '********************************************' -# # print 'HalfWidth ', HalfWidth -# # print 'Maximun ', Maximun -# # print 'eMinus1 ', eMinus1 -# # print 'Rangpos ', Rangpos -# # print 'GaussCenter ',GaussCenter -# # print 'E01 ',E01 -# # print 'N01 ',N01 -# # print 'E02 ',E02 -# # print 'N02 ',N02 -# # print 'E12 ',E12 -# # print 'N12 ',N12 -# #print 'self.dataOut.velocityX ', self.dataOut.velocityX -# # print 'Fij ', Fij -# # print 'cC ', cC -# # print 'cF ', cF -# # print 'cG ', cG -# # print 'cA ', cA -# # print 'cB ', cB -# # print 'cH ', cH -# # print 'Vx ', Vx -# # print 'Vy ', Vy -# # print 'Vmag ', Vmag -# # print 'Vang ', Vang*180/numpy.pi -# # print 'PhaseSlope ',PhaseSlope[0] -# # print 'PhaseSlope ',PhaseSlope[1] -# # print 'PhaseSlope ',PhaseSlope[2] -# # print '********************************************' -# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY) -# -# #print 'self.dataOut.velocityX', len(self.dataOut.velocityX) -# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY) -# #print 'self.dataOut.velocityV', self.dataOut.velocityV -# -# self.data_output[0]=numpy.array(self.dataOut.velocityX) -# self.data_output[1]=numpy.array(self.dataOut.velocityY) -# self.data_output[2]=numpy.array(self.dataOut.velocityV) -# -# prin= self.data_output[0][~numpy.isnan(self.data_output[0])] -# print ' ' -# print 'VmagAverage',numpy.mean(prin) -# print ' ' -# # plt.figure(5) -# # plt.subplot(211) -# # plt.plot(self.dataOut.velocityX,'yo:') -# # plt.subplot(212) -# # plt.plot(self.dataOut.velocityY,'yo:') -# -# # plt.figure(1) -# # # plt.subplot(121) -# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0') -# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1') -# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # # plt.legend() -# # plt.title('DATOS A ALTURA DE 2850 METROS') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(10) -# # # plt.subplot(121) -# # plt.plot(xFrec,ySamples[0],'b',label='Ch0') -# # plt.plot(xFrec,ySamples[1],'y',label='Ch1') -# # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # plt.legend() -# # plt.title('SELFSPECTRA EN CANALES') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # # plt.plot(xFrec,FactNorm) -# # # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(9) -# # -# # -# # plt.title('DATOS SUAVIZADOS') -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.plot(xFrec,SmoothSPC,'g') -# # -# # #plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # # -# # plt.figure(2) -# # # #plt.subplot(121) -# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra') -# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano') -# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo') -# # # #plt.plot(xFrec,phase) -# # # plt.xlabel('Suavizado, promediado KHz') -# # plt.title('SELFSPECTRA PROMEDIADO') -# # # #plt.subplot(122) -# # # #plt.plot(xSamples,zline) -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.legend() -# # # -# # # plt.figure(3) -# # # plt.subplot(311) -# # # #plt.plot(xFrec,phase[0]) -# # # plt.plot(xFrec,phase[0],'g') -# # # plt.subplot(312) -# # # plt.plot(xFrec,phase[1],'g') -# # # plt.subplot(313) -# # # plt.plot(xFrec,phase[2],'g') -# # # #plt.plot(xFrec,phase[2]) -# # # -# # # plt.figure(4) -# # # -# # # plt.plot(xSamples,coherence[0],'b') -# # # plt.plot(xSamples,coherence[1],'r') -# # # plt.plot(xSamples,coherence[2],'g') -# # plt.show() -# # # -# # # plt.clf() -# # # plt.cla() -# # # plt.close() -# -# print ' ' - - self.BlockCounter += 2 - + self.BlockCounter+=2 + else: - self.fileSelector += 1 - self.BlockCounter = 0 - print("Next File") \ No newline at end of file + self.fileSelector+=1 + self.BlockCounter=0 diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 2845d8a..322769a 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -179,9 +179,6 @@ class ParamReader(JRODataReader,ProcessingUnit): print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime)) print() -# for i in range(len(filenameList)): -# print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) - self.filenameList = filenameList self.datetimeList = datetimeList @@ -504,20 +501,11 @@ class ParamReader(JRODataReader,ProcessingUnit): def getData(self): -# if self.flagNoMoreFiles: -# self.dataOut.flagNoData = True -# print 'Process finished' -# return 0 -# if self.blockIndex==self.blocksPerFile: if not( self.__setNextFileOffline() ): self.dataOut.flagNoData = True return 0 -# if self.datablock == None: # setear esta condicion cuando no hayan datos por leers -# self.dataOut.flagNoData = True -# return 0 -# self.__readData() self.__setDataOut() self.dataOut.flagNoData = False @@ -637,7 +625,10 @@ class ParamWriter(Operation): dsDict['variable'] = self.dataList[i] #--------------------- Conditionals ------------------------ #There is no data + + if dataAux is None: + return 0 #Not array, just a number @@ -821,7 +812,7 @@ class ParamWriter(Operation): return False def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -1095,7 +1086,6 @@ class ParamWriter(Operation): return self.isConfig = True -# self.putMetadata() self.setNextFile() self.putData() diff --git a/schainpy/model/io/jroIO_spectra.py b/schainpy/model/io/jroIO_spectra.py index ea4b4af..6c523d7 100644 --- a/schainpy/model/io/jroIO_spectra.py +++ b/schainpy/model/io/jroIO_spectra.py @@ -413,9 +413,7 @@ class SpectraWriter(JRODataWriter, Operation): data_dc = None -# dataOut = None - - def __init__(self):#, **kwargs): + def __init__(self): """ Inicializador de la clase SpectraWriter para la escritura de datos de espectros. @@ -429,9 +427,7 @@ class SpectraWriter(JRODataWriter, Operation): Return: None """ - Operation.__init__(self)#, **kwargs) - - #self.isConfig = False + Operation.__init__(self) self.nTotalBlocks = 0 @@ -496,7 +492,7 @@ class SpectraWriter(JRODataWriter, Operation): def writeBlock(self): - """ + """processingHeaderObj Escribe el buffer en el file designado Affected: @@ -519,8 +515,10 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_cspc is not None: - data = numpy.zeros( self.shape_cspc_Buffer, self.dtype ) + cspc = numpy.transpose( self.data_cspc, (0,2,1) ) + #data = numpy.zeros( numpy.shape(cspc), self.dtype ) + #print 'data.shape', self.shape_cspc_Buffer if not self.processingHeaderObj.shif_fft: cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones data['real'] = cspc.real @@ -529,8 +527,9 @@ class SpectraWriter(JRODataWriter, Operation): data.tofile(self.fp) if self.data_dc is not None: - data = numpy.zeros( self.shape_dc_Buffer, self.dtype ) + dc = self.data_dc + data = numpy.zeros( numpy.shape(dc), self.dtype ) data['real'] = dc.real data['imag'] = dc.imag data = data.reshape((-1)) diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py old mode 100644 new mode 100755 index 737ed9e..5c4e2f1 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -10,11 +10,7 @@ import importlib import itertools from multiprocessing import Pool, TimeoutError from multiprocessing.pool import ThreadPool -import types -from functools import partial import time -#from sklearn.cluster import KMeans - from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters from .jroproc_base import ProcessingUnit, Operation, MPDecorator @@ -128,6 +124,7 @@ class ParametersProc(ProcessingUnit): self.dataOut.abscissaList = self.dataIn.getVelRange(1) self.dataOut.spc_noise = self.dataIn.getNoise() self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1)) + # self.dataOut.normFactor = self.dataIn.normFactor self.dataOut.pairsList = self.dataIn.pairsList self.dataOut.groupList = self.dataIn.pairsList self.dataOut.flagNoData = False @@ -136,9 +133,9 @@ class ParametersProc(ProcessingUnit): self.dataOut.ChanDist = self.dataIn.ChanDist else: self.dataOut.ChanDist = None - if hasattr(self.dataIn, 'VelRange'): #Velocities range - self.dataOut.VelRange = self.dataIn.VelRange - else: self.dataOut.VelRange = None + #if hasattr(self.dataIn, 'VelRange'): #Velocities range + # self.dataOut.VelRange = self.dataIn.VelRange + #else: self.dataOut.VelRange = None if hasattr(self.dataIn, 'RadarConst'): #Radar Constant self.dataOut.RadarConst = self.dataIn.RadarConst @@ -184,9 +181,112 @@ class ParametersProc(ProcessingUnit): def target(tups): obj, args = tups - #print 'TARGETTT', obj, args + return obj.FitGau(args) + +class SpectralFilters(Operation): + + '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR + + LimitR : It is the limit in m/s of Rainfall + LimitW : It is the limit in m/s for Winds + + Input: + + self.dataOut.data_pre : SPC and CSPC + self.dataOut.spc_range : To select wind and rainfall velocities + + Affected: + + self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind + self.dataOut.spcparam_range : Used in SpcParamPlot + self.dataOut.SPCparam : Used in PrecipitationProc + + + ''' + + def __init__(self): + Operation.__init__(self) + self.i=0 + + def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5): + + + #Limite de vientos + LimitR = PositiveLimit + LimitN = NegativeLimit + + self.spc = dataOut.data_pre[0].copy() + self.cspc = dataOut.data_pre[1].copy() + + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + VelRange = dataOut.spc_range[2] + TimeRange = dataOut.spc_range[1] + FrecRange = dataOut.spc_range[0] + + Vmax= 2*numpy.max(dataOut.spc_range[2]) + Tmax= 2*numpy.max(dataOut.spc_range[1]) + Fmax= 2*numpy.max(dataOut.spc_range[0]) + + Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()] + Breaker1R=numpy.where(VelRange == Breaker1R) + + Delta = self.Num_Bin/2 - Breaker1R[0] + + + '''Reacomodando SPCrange''' + + VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0) + + VelRange[-(self.Num_Bin/2):]+= Vmax + + FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0) + + FrecRange[-(self.Num_Bin/2):]+= Fmax + + TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0) + + TimeRange[-(self.Num_Bin/2):]+= Tmax + + ''' ------------------ ''' + + Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()] + Breaker2R=numpy.where(VelRange == Breaker2R) + + + SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1) + + SPCcut = SPCroll.copy() + for i in range(self.Num_Chn): + + SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i] + SPCcut[i,-int(Delta):,:] = dataOut.noise[i] + + SPCcut[i]=SPCcut[i]- dataOut.noise[i] + SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20 + + SPCroll[i]=SPCroll[i]-dataOut.noise[i] + SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20 + + SPC_ch1 = SPCroll + + SPC_ch2 = SPCcut + + SPCparam = (SPC_ch1, SPC_ch2, self.spc) + dataOut.SPCparam = numpy.asarray(SPCparam) + + + dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1]) + + dataOut.spcparam_range[2]=VelRange + dataOut.spcparam_range[1]=TimeRange + dataOut.spcparam_range[0]=FrecRange + + class GaussianFit(Operation): ''' @@ -198,15 +298,15 @@ class GaussianFit(Operation): self.dataOut.data_pre : SelfSpectra Output: - self.dataOut.GauSPC : SPC_ch1, SPC_ch2 + self.dataOut.SPCparam : SPC_ch1, SPC_ch2 ''' - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) self.i=0 - def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points + def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points """This routine will find a couple of generalized Gaussians to a power spectrum input: spc output: @@ -214,31 +314,12 @@ class GaussianFit(Operation): """ self.spc = dataOut.data_pre[0].copy() - - - print('SelfSpectra Shape', numpy.asarray(self.spc).shape) - - - #plt.figure(50) - #plt.subplot(121) - #plt.plot(self.spc,'k',label='spc(66)') - #plt.plot(xFrec,ySamples[1],'g',label='Ch1') - #plt.plot(xFrec,ySamples[2],'r',label='Ch2') - #plt.plot(xFrec,FitGauss,'yo:',label='fit') - #plt.legend() - #plt.title('DATOS A ALTURA DE 7500 METROS') - #plt.show() - self.Num_Hei = self.spc.shape[2] - #self.Num_Bin = len(self.spc) self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] - Vrange = dataOut.abscissaList - #print 'self.spc2', numpy.asarray(self.spc).shape - - GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei]) + GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei]) SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = numpy.NaN @@ -250,272 +331,12 @@ class GaussianFit(Operation): noise_ = dataOut.spc_noise[0].copy() - pool = Pool(processes=self.Num_Chn) args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)] objs = [self for __ in range(self.Num_Chn)] attrs = list(zip(objs, args)) gauSPC = pool.map(target, attrs) - dataOut.GauSPC = numpy.asarray(gauSPC) -# ret = [] -# for n in range(self.Num_Chn): -# self.FitGau(args[n]) -# dataOut.GauSPC = ret - - - -# for ch in range(self.Num_Chn): -# -# for ht in range(self.Num_Hei): -# #print (numpy.asarray(self.spc).shape) -# spc = numpy.asarray(self.spc)[ch,:,ht] -# -# ############################################# -# # normalizing spc and noise -# # This part differs from gg1 -# spc_norm_max = max(spc) -# spc = spc / spc_norm_max -# pnoise = pnoise / spc_norm_max -# ############################################# -# -# if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's -# fatspectra=1.0 -# else: -# fatspectra=0.5 -# -# wnoise = noise_ / spc_norm_max -# #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise -# #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used -# #if wnoise>1.1*pnoise: # to be tested later -# # wnoise=pnoise -# noisebl=wnoise*0.9; noisebh=wnoise*1.1 -# spc=spc-wnoise -# -# minx=numpy.argmin(spc) -# spcs=numpy.roll(spc,-minx) -# cum=numpy.cumsum(spcs) -# tot_noise=wnoise * self.Num_Bin #64; -# #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' -# #snr=tot_signal/tot_noise -# #snr=cum[-1]/tot_noise -# -# #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise -# -# snr = sum(spcs)/tot_noise -# snrdB=10.*numpy.log10(snr) -# -# #if snrdB < -9 : -# # snrdB = numpy.NaN -# # continue -# -# #print 'snr',snrdB # , sum(spcs) , tot_noise -# -# -# #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: -# # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None -# -# cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region -# cumlo=cummax*epsi; -# cumhi=cummax*(1-epsi) -# powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-9: # when SNR is strong pick the peak with least shift (LOS velocity) error -# if oneG: -# choice=0 -# else: -# w1=lsq2[0][1]; w2=lsq2[0][5] -# a1=lsq2[0][2]; a2=lsq2[0][6] -# p1=lsq2[0][3]; p2=lsq2[0][7] -# s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; -# gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling -# -# if gp1>gp2: -# if a1>0.7*a2: -# choice=1 -# else: -# choice=2 -# elif gp2>gp1: -# if a2>0.7*a1: -# choice=2 -# else: -# choice=1 -# else: -# choice=numpy.argmax([a1,a2])+1 -# #else: -# #choice=argmin([std2a,std2b])+1 -# -# else: # with low SNR go to the most energetic peak -# choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) -# -# #print 'choice',choice -# -# if choice==0: # pick the single gaussian fit -# Amplitude0=lsq1[0][2] -# shift0=lsq1[0][0] -# width0=lsq1[0][1] -# p0=lsq1[0][3] -# Amplitude1=0. -# shift1=0. -# width1=0. -# p1=0. -# noise=lsq1[0][4] -# elif choice==1: # take the first one of the 2 gaussians fitted -# Amplitude0 = lsq2[0][2] -# shift0 = lsq2[0][0] -# width0 = lsq2[0][1] -# p0 = lsq2[0][3] -# Amplitude1 = lsq2[0][6] # This is 0 in gg1 -# shift1 = lsq2[0][4] # This is 0 in gg1 -# width1 = lsq2[0][5] # This is 0 in gg1 -# p1 = lsq2[0][7] # This is 0 in gg1 -# noise = lsq2[0][8] -# else: # the second one -# Amplitude0 = lsq2[0][6] -# shift0 = lsq2[0][4] -# width0 = lsq2[0][5] -# p0 = lsq2[0][7] -# Amplitude1 = lsq2[0][2] # This is 0 in gg1 -# shift1 = lsq2[0][0] # This is 0 in gg1 -# width1 = lsq2[0][1] # This is 0 in gg1 -# p1 = lsq2[0][3] # This is 0 in gg1 -# noise = lsq2[0][8] -# -# #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) -# SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 -# SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 -# #print 'SPC_ch1.shape',SPC_ch1.shape -# #print 'SPC_ch2.shape',SPC_ch2.shape -# #dataOut.data_param = SPC_ch1 -# GauSPC[0] = SPC_ch1 -# GauSPC[1] = SPC_ch2 - -# #plt.gcf().clear() -# plt.figure(50+self.i) -# self.i=self.i+1 -# #plt.subplot(121) -# plt.plot(self.spc,'k')#,label='spc(66)') -# plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1') -# #plt.plot(SPC_ch2,'r')#,label='gg2') -# #plt.plot(xFrec,ySamples[1],'g',label='Ch1') -# #plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# #plt.plot(xFrec,FitGauss,'yo:',label='fit') -# plt.legend() -# plt.title('DATOS A ALTURA DE 7500 METROS') -# plt.show() -# print 'shift0', shift0 -# print 'Amplitude0', Amplitude0 -# print 'width0', width0 -# print 'p0', p0 -# print '========================' -# print 'shift1', shift1 -# print 'Amplitude1', Amplitude1 -# print 'width1', width1 -# print 'p1', p1 -# print 'noise', noise -# print 's_noise', wnoise - - print('========================================================') - print('total_time: ', time.time()-start_time) - - # re-normalizing spc and noise - # This part differs from gg1 - - + dataOut.SPCparam = numpy.asarray(SPCparam) ''' Parameters: 1. Amplitude @@ -524,16 +345,11 @@ class GaussianFit(Operation): 4. Power ''' - - ############################################################################### def FitGau(self, X): Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X - #print 'VARSSSS', ch, pnoise, noise, num_intg - - #print 'HEIGHTS', self.Num_Hei - - GauSPC = [] + + SPCparam = [] SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch1[:] = 0#numpy.NaN @@ -542,10 +358,6 @@ class GaussianFit(Operation): for ht in range(self.Num_Hei): - #print (numpy.asarray(self.spc).shape) - - #print 'TTTTT', ch , ht - #print self.spc.shape spc = numpy.asarray(self.spc)[ch,:,ht] @@ -554,27 +366,26 @@ class GaussianFit(Operation): # normalizing spc and noise # This part differs from gg1 spc_norm_max = max(spc) - spc = spc / spc_norm_max - pnoise = pnoise / spc_norm_max + #spc = spc / spc_norm_max + pnoise = pnoise #/ spc_norm_max ############################################# - + fatspectra=1.0 - wnoise = noise_ / spc_norm_max + wnoise = noise_ #/ spc_norm_max #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used #if wnoise>1.1*pnoise: # to be tested later # wnoise=pnoise - noisebl=wnoise*0.9; noisebh=wnoise*1.1 + noisebl=wnoise*0.9; + noisebh=wnoise*1.1 spc=spc-wnoise - # print 'wnoise', noise_[0], spc_norm_max, wnoise + minx=numpy.argmin(spc) + #spcs=spc.copy() spcs=numpy.roll(spc,-minx) cum=numpy.cumsum(spcs) tot_noise=wnoise * self.Num_Bin #64; - #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise - #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' - #snr=tot_signal/tot_noise - #snr=cum[-1]/tot_noise + snr = sum(spcs)/tot_noise snrdB=10.*numpy.log10(snr) @@ -582,16 +393,15 @@ class GaussianFit(Operation): snr = numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN SPC_ch1[:,ht] = 0#numpy.NaN - GauSPC = (SPC_ch1,SPC_ch2) + SPCparam = (SPC_ch1,SPC_ch2) continue - #print 'snr',snrdB #, sum(spcs) , tot_noise - #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None - cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region + cummax=max(cum); + epsi=0.08*fatspectra # cumsum to narrow down the energy region cumlo=cummax*epsi; cumhi=cummax*(1-epsi) powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum-6: # when SNR is strong pick the peak with least shift (LOS velocity) error + + if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error if oneG: choice=0 else: @@ -690,7 +483,7 @@ class GaussianFit(Operation): s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2; gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling - + if gp1>gp2: if a1>0.7*a2: choice=1 @@ -710,13 +503,15 @@ class GaussianFit(Operation): choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]]) - shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) - shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) + shift0=lsq2[0][0]; + vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0]) + shift1=lsq2[0][4]; + vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0]) - max_vel = 20 + max_vel = 1.0 #first peak will be 0, second peak will be 1 - if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range + if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range shift0=lsq2[0][0] width0=lsq2[0][1] Amplitude0=lsq2[0][2] @@ -739,120 +534,18 @@ class GaussianFit(Operation): p0=lsq2[0][7] noise=lsq2[0][8] - if Amplitude0<0.1: # in case the peak is noise - shift0,width0,Amplitude0,p0 = 4*[numpy.NaN] - if Amplitude1<0.1: - shift1,width1,Amplitude1,p1 = 4*[numpy.NaN] - - -# if choice==0: # pick the single gaussian fit -# Amplitude0=lsq1[0][2] -# shift0=lsq1[0][0] -# width0=lsq1[0][1] -# p0=lsq1[0][3] -# Amplitude1=0. -# shift1=0. -# width1=0. -# p1=0. -# noise=lsq1[0][4] -# elif choice==1: # take the first one of the 2 gaussians fitted -# Amplitude0 = lsq2[0][2] -# shift0 = lsq2[0][0] -# width0 = lsq2[0][1] -# p0 = lsq2[0][3] -# Amplitude1 = lsq2[0][6] # This is 0 in gg1 -# shift1 = lsq2[0][4] # This is 0 in gg1 -# width1 = lsq2[0][5] # This is 0 in gg1 -# p1 = lsq2[0][7] # This is 0 in gg1 -# noise = lsq2[0][8] -# else: # the second one -# Amplitude0 = lsq2[0][6] -# shift0 = lsq2[0][4] -# width0 = lsq2[0][5] -# p0 = lsq2[0][7] -# Amplitude1 = lsq2[0][2] # This is 0 in gg1 -# shift1 = lsq2[0][0] # This is 0 in gg1 -# width1 = lsq2[0][1] # This is 0 in gg1 -# p1 = lsq2[0][3] # This is 0 in gg1 -# noise = lsq2[0][8] - - #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) + if Amplitude0<0.05: # in case the peak is noise + shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN] + if Amplitude1<0.05: + shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] + + SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 - #print 'SPC_ch1.shape',SPC_ch1.shape - #print 'SPC_ch2.shape',SPC_ch2.shape - #dataOut.data_param = SPC_ch1 - GauSPC = (SPC_ch1,SPC_ch2) - #GauSPC[1] = SPC_ch2 - -# print 'shift0', shift0 -# print 'Amplitude0', Amplitude0 -# print 'width0', width0 -# print 'p0', p0 -# print '========================' -# print 'shift1', shift1 -# print 'Amplitude1', Amplitude1 -# print 'width1', width1 -# print 'p1', p1 -# print 'noise', noise -# print 's_noise', wnoise + SPCparam = (SPC_ch1,SPC_ch2) - return GauSPC - - - def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination. - y_model=self.y_model1(x,state) - s0,w0,a0,p0,n=state - e0=((x-s0)/w0)**2; - - e0u=((x-s0-self.Num_Bin)/w0)**2; - - e0d=((x-s0+self.Num_Bin)/w0)**2 - m0=numpy.exp(-0.5*e0**(p0/2.)); - m0u=numpy.exp(-0.5*e0u**(p0/2.)); - m0d=numpy.exp(-0.5*e0d**(p0/2.)) - JA=m0+m0u+m0d - JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) - - JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) - - JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 - jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model]) - return jack1.T - - def y_jacobian2(self,x,state): - y_model=self.y_model2(x,state) - s0,w0,a0,p0,s1,w1,a1,p1,n=state - e0=((x-s0)/w0)**2; - - e0u=((x-s0- self.Num_Bin )/w0)**2; - - e0d=((x-s0+ self.Num_Bin )/w0)**2 - e1=((x-s1)/w1)**2; - e1u=((x-s1- self.Num_Bin )/w1)**2; - - e1d=((x-s1+ self.Num_Bin )/w1)**2 - m0=numpy.exp(-0.5*e0**(p0/2.)); - m0u=numpy.exp(-0.5*e0u**(p0/2.)); - m0d=numpy.exp(-0.5*e0d**(p0/2.)) - m1=numpy.exp(-0.5*e1**(p1/2.)); - m1u=numpy.exp(-0.5*e1u**(p1/2.)); - m1d=numpy.exp(-0.5*e1d**(p1/2.)) - JA=m0+m0u+m0d - JA1=m1+m1u+m1d - JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d) - JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d) - - JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0) - - JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1) - - JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2 - - JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2 - jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model]) - return jack2.T + return GauSPC def y_model1(self,x,state): shift0,width0,amplitude0,power0,noise=state @@ -884,6 +577,7 @@ class GaussianFit(Operation): def misfit2(self,state,y_data,x,num_intg): return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) + class PrecipitationProc(Operation): @@ -900,24 +594,61 @@ class PrecipitationProc(Operation): Parameters affected: ''' - - def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None, - tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None): + def __init__(self): + Operation.__init__(self) + self.i=0 + + + def gaus(self,xSamples,Amp,Mu,Sigma): + return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) + + + + def Moments(self, ySamples, xSamples): + Pot = numpy.nansum( ySamples ) # Potencia, momento 0 + yNorm = ySamples / Pot + + Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral + + return numpy.array([Pot, Vr, Desv]) + + def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118, + tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350): - self.spc = dataOut.data_pre[0].copy() - self.Num_Hei = self.spc.shape[2] - self.Num_Bin = self.spc.shape[1] - self.Num_Chn = self.spc.shape[0] - Velrange = dataOut.abscissaList + Velrange = dataOut.spcparam_range[2] + FrecRange = dataOut.spcparam_range[0] + + dV= Velrange[1]-Velrange[0] + dF= FrecRange[1]-FrecRange[0] if radar == "MIRA35C" : + self.spc = dataOut.data_pre[0].copy() + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] Ze = self.dBZeMODE2(dataOut) else: + self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() # + + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + self.spc[:,:,0:7]= numpy.NaN + + """##########################################""" + + self.Num_Hei = self.spc.shape[2] + self.Num_Bin = self.spc.shape[1] + self.Num_Chn = self.spc.shape[0] + + ''' Se obtiene la constante del RADAR ''' + self.Pt = Pt self.Gt = Gt self.Gr = Gr @@ -927,48 +658,101 @@ class PrecipitationProc(Operation): self.ThetaT = ThetaT self.ThetaR = ThetaR - RadarConstant = GetRadarConstant() - SPCmean = numpy.mean(self.spc,0) - ETA = numpy.zeros(self.Num_Hei) - Pr = numpy.sum(SPCmean,0) + Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) + Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR) + RadarConstant = 5e-26 * Numerator / Denominator # - #for R in range(self.Num_Hei): - # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) - - D_range = numpy.zeros(self.Num_Hei) - EqSec = numpy.zeros(self.Num_Hei) + ''' ============================= ''' + + self.spc[0] = (self.spc[0]-dataOut.noise[0]) + self.spc[1] = (self.spc[1]-dataOut.noise[1]) + self.spc[2] = (self.spc[2]-dataOut.noise[2]) + + self.spc[ numpy.where(self.spc < 0)] = 0 + + SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise)) + SPCmean[ numpy.where(SPCmean < 0)] = 0 + + ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei]) + ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei]) + ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei]) + + Pr = SPCmean[:,:] + + VelMeteoro = numpy.mean(SPCmean,axis=0) + + D_range = numpy.zeros([self.Num_Bin,self.Num_Hei]) + SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei]) + N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei]) + V_mean = numpy.zeros(self.Num_Hei) del_V = numpy.zeros(self.Num_Hei) + Z = numpy.zeros(self.Num_Hei) + Ze = numpy.zeros(self.Num_Hei) + RR = numpy.zeros(self.Num_Hei) + + Range = dataOut.heightList*1000. for R in range(self.Num_Hei): - ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA) - h = R + Altitude #Range from ground to radar pulse altitude + h = Range[R] + Altitude #Range from ground to radar pulse altitude del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity - D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity - SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma) + D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3 + + '''NOTA: ETA(n) dn = ETA(f) df + + dn = 1 Diferencial de muestreo + df = ETA(n) / ETA(f) + + ''' + + ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA) + + ETAv[:,R]=ETAn[:,R]/dV + + ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R]) + + SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma) + + N_dist[:,R] = ETAn[:,R] / SIGMA[:,R] + + DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin]) + + try: + popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments) + except: + popt01=numpy.zeros(3) + popt01[1]= DMoments[1] + + if popt01[1]<0 or popt01[1]>20: + popt01[1]=numpy.NaN + + + V_mean[R]=popt01[1] + + Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18 + + RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate + + Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km) - N_dist[R] = ETA[R] / SIGMA[R] - - Ze = (ETA * Lambda**4) / (numpy.pi * Km) - Z = numpy.sum( N_dist * D_range**6 ) - RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate - RR = (Ze/200)**(1/1.6) + RR2 = (Z/200)**(1/1.6) dBRR = 10*numpy.log10(RR) + dBRR2 = 10*numpy.log10(RR2) dBZe = 10*numpy.log10(Ze) - dataOut.data_output = Ze - dataOut.data_param = numpy.ones([2,self.Num_Hei]) - dataOut.channelList = [0,1] - print('channelList', dataOut.channelList) - dataOut.data_param[0]=dBZe - dataOut.data_param[1]=dBRR - print('RR SHAPE', dBRR.shape) - print('Ze SHAPE', dBZe.shape) - print('dataOut.data_param SHAPE', dataOut.data_param.shape) - + dBZ = 10*numpy.log10(Z) + + dataOut.data_output = RR[8] + dataOut.data_param = numpy.ones([3,self.Num_Hei]) + dataOut.channelList = [0,1,2] + + dataOut.data_param[0]=dBZ + dataOut.data_param[1]=V_mean + dataOut.data_param[2]=RR + def dBZeMODE2(self, dataOut): # Processing for MIRA35C @@ -983,7 +767,7 @@ class PrecipitationProc(Operation): data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN ETA = numpy.sum(SNR,1) - print('ETA' , ETA) + ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN) Ze = numpy.ones([self.Num_Chn, self.Num_Hei] ) @@ -995,26 +779,27 @@ class PrecipitationProc(Operation): return Ze - def GetRadarConstant(self): - - """ - Constants: - - Pt: Transmission Power dB - Gt: Transmission Gain dB - Gr: Reception Gain dB - Lambda: Wavelenght m - aL: Attenuation loses dB - tauW: Width of transmission pulse s - ThetaT: Transmission antenna bean angle rad - ThetaR: Reception antenna beam angle rad - - """ - Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) - Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) - RadarConstant = Numerator / Denominator - - return RadarConstant +# def GetRadarConstant(self): +# +# """ +# Constants: +# +# Pt: Transmission Power dB 5kW 5000 +# Gt: Transmission Gain dB 24.7 dB 295.1209 +# Gr: Reception Gain dB 18.5 dB 70.7945 +# Lambda: Wavelenght m 0.6741 m 0.6741 +# aL: Attenuation loses dB 4dB 2.5118 +# tauW: Width of transmission pulse s 4us 4e-6 +# ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317 +# ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087 +# +# """ +# +# Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) +# Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR) +# RadarConstant = Numerator / Denominator +# +# return RadarConstant @@ -1037,10 +822,20 @@ class FullSpectralAnalysis(Operation): Parameters affected: Winds, height range, SNR """ - def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7): + def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7): + + self.indice=int(numpy.random.rand()*1000) spc = dataOut.data_pre[0].copy() - cspc = dataOut.data_pre[1].copy() + cspc = dataOut.data_pre[1] + + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + SNRspc = spc.copy() + SNRspc[:,:,0:7]= numpy.NaN + + """##########################################""" + nChannel = spc.shape[0] nProfiles = spc.shape[1] @@ -1050,14 +845,9 @@ class FullSpectralAnalysis(Operation): if dataOut.ChanDist is not None : ChanDist = dataOut.ChanDist else: - ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]]) - - #print 'ChanDist', ChanDist + ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]]) - if dataOut.VelRange is not None: - VelRange= dataOut.VelRange - else: - VelRange= dataOut.abscissaList + FrecRange = dataOut.spc_range[0] ySamples=numpy.ones([nChannel,nProfiles]) phase=numpy.ones([nChannel,nProfiles]) @@ -1065,113 +855,108 @@ class FullSpectralAnalysis(Operation): coherence=numpy.ones([nChannel,nProfiles]) PhaseSlope=numpy.ones(nChannel) PhaseInter=numpy.ones(nChannel) - dataSNR = dataOut.data_SNR - - + data_SNR=numpy.zeros([nProfiles]) data = dataOut.data_pre noise = dataOut.noise - print('noise',noise) - #SNRdB = 10*numpy.log10(dataOut.data_SNR) - FirstMoment = numpy.average(dataOut.data_param[:,1,:],0) - #SNRdBMean = [] - + dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0] - #for j in range(nHeights): - # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]])) - # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]])) - - data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN + dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20 + + + data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN velocityX=[] velocityY=[] velocityV=[] + PhaseLine=[] - dbSNR = 10*numpy.log10(dataSNR) + dbSNR = 10*numpy.log10(dataOut.data_SNR) dbSNR = numpy.average(dbSNR,0) + for Height in range(nHeights): - - [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit) + + [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit) + PhaseLine = numpy.append(PhaseLine, PhaseSlope) if abs(Vzon)<100. and abs(Vzon)> 0.: velocityX=numpy.append(velocityX, Vzon)#Vmag else: - print('Vzon',Vzon) velocityX=numpy.append(velocityX, numpy.NaN) if abs(Vmer)<100. and abs(Vmer) > 0.: - velocityY=numpy.append(velocityY, Vmer)#Vang + velocityY=numpy.append(velocityY, -Vmer)#Vang else: - print('Vmer',Vmer) velocityY=numpy.append(velocityY, numpy.NaN) if dbSNR[Height] > SNRlimit: - velocityV=numpy.append(velocityV, FirstMoment[Height]) + velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height]) else: velocityV=numpy.append(velocityV, numpy.NaN) - #FirstMoment[Height]= numpy.NaN -# if SNRdBMean[Height] <12: -# FirstMoment[Height] = numpy.NaN -# velocityX[Height] = numpy.NaN -# velocityY[Height] = numpy.NaN - - - data_output[0]=numpy.array(velocityX) - data_output[1]=numpy.array(velocityY) - data_output[2]=-velocityV#FirstMoment - - print(' ') - #print 'FirstMoment' - #print FirstMoment - print('velocityX',data_output[0]) - print(' ') - print('velocityY',data_output[1]) - #print numpy.array(velocityY) - print(' ') - #print 'SNR' - #print 10*numpy.log10(dataOut.data_SNR) - #print numpy.shape(10*numpy.log10(dataOut.data_SNR)) - print(' ') + + + '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR''' + data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1) + data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1) + data_output[2] = velocityV#FirstMoment + + xFrec=FrecRange[0:spc.shape[1]] dataOut.data_output=data_output + return def moving_average(self,x, N=2): return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):] - def gaus(self,xSamples,a,x0,sigma): - return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2)) + def gaus(self,xSamples,Amp,Mu,Sigma): + return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) + + - def Find(self,x,value): - for index in range(len(x)): - if x[index]==value: - return index + def Moments(self, ySamples, xSamples): + Pot = numpy.nansum( ySamples ) # Potencia, momento 0 + yNorm = ySamples / Pot + Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento + Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento + Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral + + return numpy.array([Pot, Vr, Desv]) - def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit): + def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit): + + ySamples=numpy.ones([spc.shape[0],spc.shape[1]]) phase=numpy.ones([spc.shape[0],spc.shape[1]]) CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) coherence=numpy.ones([spc.shape[0],spc.shape[1]]) - PhaseSlope=numpy.ones(spc.shape[0]) + PhaseSlope=numpy.zeros(spc.shape[0]) PhaseInter=numpy.ones(spc.shape[0]) - xFrec=VelRange + xFrec=AbbsisaRange[0][0:spc.shape[1]] + xVel =AbbsisaRange[2][0:spc.shape[1]] + Vv=numpy.empty(spc.shape[2])*0 + SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]# + + SPCmoments = self.Moments(SPCav[:,Height], xVel ) + CSPCmoments = [] + cspcNoise = numpy.empty(3) '''Getting Eij and Nij''' - E01=ChanDist[0][0] - N01=ChanDist[0][1] + Xi01=ChanDist[0][0] + Eta01=ChanDist[0][1] - E02=ChanDist[1][0] - N02=ChanDist[1][1] + Xi02=ChanDist[1][0] + Eta02=ChanDist[1][1] - E12=ChanDist[2][0] - N12=ChanDist[2][1] + Xi12=ChanDist[2][0] + Eta12=ChanDist[2][1] z = spc.copy() z = numpy.where(numpy.isfinite(z), z, numpy.NAN) @@ -1179,176 +964,197 @@ class FullSpectralAnalysis(Operation): for i in range(spc.shape[0]): '''****** Line of Data SPC ******''' - zline=z[i,:,Height] + zline=z[i,:,Height].copy() - noise[i] # Se resta ruido '''****** SPC is normalized ******''' - FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy()) - FactNorm= FactNorm/numpy.sum(FactNorm) + SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido + FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado - SmoothSPC=self.moving_average(FactNorm,N=3) - - xSamples = ar(list(range(len(SmoothSPC)))) - ySamples[i] = SmoothSPC - - #dbSNR=10*numpy.log10(dataSNR) - print(' ') - print(' ') - print(' ') - - #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120] - print('SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20]) - print('noise',noise) - print('zline',zline.shape, zline[0:20]) - print('FactNorm',FactNorm.shape, FactNorm[0:20]) - print('FactNorm suma', numpy.sum(FactNorm)) + xSamples = xFrec # Se toma el rango de frecuncias + ySamples[i] = FactNorm # Se toman los valores de SPC normalizado for i in range(spc.shape[0]): '''****** Line of Data CSPC ******''' - cspcLine=cspc[i,:,Height].copy() + cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido + SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido + cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado - '''****** CSPC is normalized ******''' + '''****** CSPC is normalized with respect to Briggs and Vincent ******''' chan_index0 = pairsList[i][0] chan_index1 = pairsList[i][1] - CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) # - CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor) + CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2 + CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor) CSPCSamples[i] = CSPCNorm + coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor) - coherence[i]= self.moving_average(coherence[i],N=2) + #coherence[i]= self.moving_average(coherence[i],N=1) phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi - print('cspcLine', cspcLine.shape, cspcLine[0:20]) - print('CSPCFactor', CSPCFactor)#, CSPCFactor[0:20] - print(numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i]) - print('CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20]) - print('CSPCNorm suma', numpy.sum(CSPCNorm)) - print('CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20]) + CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples), + self.Moments(numpy.abs(CSPCSamples[1]), xSamples), + self.Moments(numpy.abs(CSPCSamples[2]), xSamples)]) - '''****** Getting fij width ******''' + + popt=[1e-10,0,1e-10] + popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10] + FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0 + + CSPCMask01 = numpy.abs(CSPCSamples[0]) + CSPCMask02 = numpy.abs(CSPCSamples[1]) + CSPCMask12 = numpy.abs(CSPCSamples[2]) + + mask01 = ~numpy.isnan(CSPCMask01) + mask02 = ~numpy.isnan(CSPCMask02) + mask12 = ~numpy.isnan(CSPCMask12) - yMean=[] - yMean2=[] + #mask = ~numpy.isnan(CSPCMask01) + CSPCMask01 = CSPCMask01[mask01] + CSPCMask02 = CSPCMask02[mask02] + CSPCMask12 = CSPCMask12[mask12] + #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01) - for j in range(len(ySamples[1])): - yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]])) - '''******* Getting fitting Gaussian ******''' - meanGauss=sum(xSamples*yMean) / len(xSamples) - sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) - print('****************************') - print('len(xSamples): ',len(xSamples)) - print('yMean: ', yMean.shape, yMean[0:20]) - print('ySamples', ySamples.shape, ySamples[0,0:20]) - print('xSamples: ',xSamples.shape, xSamples[0:20]) + '''***Fit Gauss CSPC01***''' + if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 : + try: + popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0]) + popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1]) + popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2]) + FitGauss01 = self.gaus(xSamples,*popt01) + FitGauss02 = self.gaus(xSamples,*popt02) + FitGauss12 = self.gaus(xSamples,*popt12) + except: + FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0])) + FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1])) + FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2])) + + + CSPCopt = numpy.vstack([popt01,popt02,popt12]) + + '''****** Getting fij width ******''' + + yMean = numpy.average(ySamples, axis=0) # ySamples[0] + + '''******* Getting fitting Gaussian *******''' + meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia) + sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia) - print('meanGauss',meanGauss) - print('sigma',sigma) + yMoments = self.Moments(yMean, xSamples) - #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001): - if dbSNR > SNRlimit : - try: - popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) + if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001: + try: + popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments) + FitGauss=self.gaus(xSamples,*popt) - if numpy.amax(popt)>numpy.amax(yMean)*0.3: - FitGauss=self.gaus(xSamples,*popt) - - else: - FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - print('Verificador: Dentro', Height) except :#RuntimeError: FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - + else: FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) - Maximun=numpy.amax(yMean) - eMinus1=Maximun*numpy.exp(-1)#*0.8 - HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) - HalfWidth= xFrec[HWpos] - GCpos=self.Find(FitGauss, numpy.amax(FitGauss)) - Vpos=self.Find(FactNorm, numpy.amax(FactNorm)) - - #Vpos=FirstMoment[] '''****** Getting Fij ******''' + Fijcspc = CSPCopt[:,2]/2*3 + + + GaussCenter = popt[1] #xFrec[GCpos] + #Punto en Eje X de la Gaussiana donde se encuentra el centro + ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()] + PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0] + + #Punto e^-1 hubicado en la Gaussiana + PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1) + FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss" + PointFij = numpy.where(FitGauss==FijClosest)[0][0] - GaussCenter=xFrec[GCpos] - if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): - Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 + if xSamples[PointFij] > xSamples[PointGauCenter]: + Fij = xSamples[PointFij] - xSamples[PointGauCenter] + else: - Fij=abs(GaussCenter-HalfWidth)+0.0000001 + Fij = xSamples[PointGauCenter] - xSamples[PointFij] + + + '''****** Taking frequency ranges from SPCs ******''' - '''****** Getting Frecuency range of significant data ******''' - Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) + #GaussCenter = popt[1] #Primer momento 01 + GauWidth = popt[2] *3/2 #Ancho de banda de Gau01 + Range = numpy.empty(2) + Range[0] = GaussCenter - GauWidth + Range[1] = GaussCenter + GauWidth + #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) + ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()] + ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()] - if Rangpos5 and len(FrecRange)5 and len(FrecRange)4: + Vver=popt[1] + else: + Vver=numpy.NaN + FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12]) + + + return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC + class SpectralMoments(Operation): ''' @@ -1384,7 +1195,7 @@ class SpectralMoments(Operation): self.dataOut.noise : Noise level per channel Affected: - self.dataOut.data_param : Parameters per channel + self.dataOut.moments : Parameters per channel self.dataOut.data_SNR : SNR per channel ''' @@ -1401,7 +1212,7 @@ class SpectralMoments(Operation): for ind in range(nChannel): data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] ) - dataOut.data_param = data_param[:,1:,:] + dataOut.moments = data_param[:,1:,:] dataOut.data_SNR = data_param[:,0] dataOut.data_DOP = data_param[:,1] dataOut.data_MEAN = data_param[:,2] @@ -1431,6 +1242,8 @@ class SpectralMoments(Operation): vec_fd = numpy.zeros(oldspec.shape[1]) vec_w = numpy.zeros(oldspec.shape[1]) vec_snr = numpy.zeros(oldspec.shape[1]) + + oldspec = numpy.ma.masked_invalid(oldspec) for ind in range(oldspec.shape[1]): @@ -1469,7 +1282,7 @@ class SpectralMoments(Operation): fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) snr = (spec2.mean()-n0)/n0 - + if (snr < 1.e-20) : snr = 1.e-20 @@ -1477,7 +1290,7 @@ class SpectralMoments(Operation): vec_fd[ind] = fd vec_w[ind] = w vec_snr[ind] = snr - + moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) return moments @@ -1675,7 +1488,6 @@ class SpectralFitting(Operation): dataCross = dataCross**2/K for h in range(nHeights): -# print self.dataOut.heightList[h] #Input d = data[:,h] @@ -1734,7 +1546,7 @@ class SpectralFitting(Operation): fm = self.dataOut.library.modelFunction(p, constants) fmp=numpy.dot(LT,fm) - + return dp-fmp def __getSNR(self, z, noise): @@ -1768,8 +1580,8 @@ class WindProfiler(Operation): n = None - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) def __calculateCosDir(self, elev, azim): zen = (90 - elev)*numpy.pi/180 @@ -2071,12 +1883,9 @@ class WindProfiler(Operation): Parameters affected: Winds ''' -# print arrayMeteor.shape #Settings nInt = (heightMax - heightMin)/2 -# print nInt nInt = int(nInt) -# print nInt winds = numpy.zeros((2,nInt))*numpy.nan #Filter errors @@ -2475,8 +2284,8 @@ class WindProfiler(Operation): class EWDriftsEstimation(Operation): - def __init__(self, **kwargs): - Operation.__init__(self, **kwargs) + def __init__(self): + Operation.__init__(self) def __correctValues(self, heiRang, phi, velRadial, SNR): listPhi = phi.tolist() diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index db71e0d..5c68c91 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -159,9 +159,7 @@ class SpectraProc(ProcessingUnit): dtype='complex') if self.dataIn.flagDataAsBlock: - # data dimension: [nChannels, nProfiles, nSamples] nVoltProfiles = self.dataIn.data.shape[1] - # nVoltProfiles = self.dataIn.nProfiles if nVoltProfiles == nProfiles: self.buffer = self.dataIn.data.copy() @@ -299,7 +297,57 @@ class SpectraProc(ProcessingUnit): self.__selectPairsByChannel(self.dataOut.channelList) return 1 + + + def selectFFTs(self, minFFT, maxFFT ): + """ + Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango + minFFT<= FFT <= maxFFT + """ + + if (minFFT > maxFFT): + raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT)) + + if (minFFT < self.dataOut.getFreqRange()[0]): + minFFT = self.dataOut.getFreqRange()[0] + + if (maxFFT > self.dataOut.getFreqRange()[-1]): + maxFFT = self.dataOut.getFreqRange()[-1] + + minIndex = 0 + maxIndex = 0 + FFTs = self.dataOut.getFreqRange() + + inda = numpy.where(FFTs >= minFFT) + indb = numpy.where(FFTs <= maxFFT) + + try: + minIndex = inda[0][0] + except: + minIndex = 0 + + try: + maxIndex = indb[0][-1] + except: + maxIndex = len(FFTs) + + self.selectFFTsByIndex(minIndex, maxIndex) + return 1 + + + def setH0(self, h0, deltaHeight = None): + + if not deltaHeight: + deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0] + + nHeights = self.dataOut.nHeights + + newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight + + self.dataOut.heightList = newHeiRange + + def selectHeights(self, minHei, maxHei): """ Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango @@ -316,9 +364,9 @@ class SpectraProc(ProcessingUnit): 1 si el metodo se ejecuto con exito caso contrario devuelve 0 """ + if (minHei > maxHei): - raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % ( - minHei, maxHei)) + raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)) if (minHei < self.dataOut.heightList[0]): minHei = self.dataOut.heightList[0] @@ -344,6 +392,7 @@ class SpectraProc(ProcessingUnit): maxIndex = len(heights) self.selectHeightsByIndex(minIndex, maxIndex) + return 1 @@ -389,6 +438,40 @@ class SpectraProc(ProcessingUnit): return 1 + def selectFFTsByIndex(self, minIndex, maxIndex): + """ + + """ + + if (minIndex < 0) or (minIndex > maxIndex): + raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex)) + + if (maxIndex >= self.dataOut.nProfiles): + maxIndex = self.dataOut.nProfiles-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.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1]) + self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1] + self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1] + + return 1 + + + def selectHeightsByIndex(self, minIndex, maxIndex): """ Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango @@ -494,7 +577,32 @@ class SpectraProc(ProcessingUnit): return 1 - def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None): + def removeInterference2(self): + + cspc = self.dataOut.data_cspc + spc = self.dataOut.data_spc + Heights = numpy.arange(cspc.shape[2]) + realCspc = numpy.abs(cspc) + + for i in range(cspc.shape[0]): + LinePower= numpy.sum(realCspc[i], axis=0) + Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)] + SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ] + InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 ) + InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)] + InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)] + + + InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) ) + #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax])) + if len(InterferenceRange)