schain Commit - r1385:83c0b1494611 · Jicamarca Repository

Adicionada operación Clean Rayleigh a Spectra Proc

joabAM -

r1385:83c0b1494611

parent child

Context file:

r1385:83c0b1494611

Collapse all files

schainpy/model/graphics/jroplot_heispectra.py +8 -7

             import numpy
             from schainpy.model.graphics.jroplot_base import Plot, plt
+            import matplotlib.pyplot as plt
             class SpectraHeisPlot(Plot):
                     meta = {}
                     spc = 10*numpy.log10(dataOut.data_spc / dataOut.normFactor)
                     data['spc_heis'] = spc
                     return data, meta
                 def plot(self):
                     c = 3E8
-                    deltaHeight = self.data.yrange[1] - self.data.yrange[0]
+                    deltaHeight = self.data.yrange[1] - self.data.yrange[0] # yrange = heightList
                     x = numpy.arange(-1*len(self.data.yrange)/2., len(self.data.yrange)/2.)*(c/(2*deltaHeight*len(self.data.yrange)*1000))
                     self.y = self.data[-1]['spc_heis']
                     self.titles = []
+                    Maintitle = "Range from %d km to %d km" %(int(self.data.yrange[0]),int(self.data.yrange[-1]))
                     for n, ax in enumerate(self.axes):
                         ychannel = self.y[n,:]
                         if ax.firsttime:
                             ax.plt.set_data(x, ychannel)
                         self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
+                    plt.suptitle(Maintitle)
             class RTIHeisPlot(Plot):
                     spc = dataOut.data_spc / dataOut.normFactor
                     spc = 10*numpy.log10(numpy.average(spc, axis=1))
                     data['rti_heis'] = spc
                     return data, meta
                 def plot(self):

schainpy/model/graphics/jroplot_voltage.py +7 -7

                 def setup(self):
                     self.xaxis = 'Range (Km)'
-                    self.ncols = 1
+                    self.nplots = len(self.data.channels)
-                    self.nrows = 1
+                    self.nrows = int(numpy.ceil(self.nplots/2))
-                    self.nplots = 1
+                    self.ncols = int(numpy.ceil(self.nplots/self.nrows))
                     self.ylabel = 'Intensity [dB]'
-                    self.titles = ['Scope']
+                    self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
                     self.colorbar = False
                     self.width = 6
                     self.height = 4
                     yreal = y[channelIndexList,:].real
                     yimag = y[channelIndexList,:].imag
-                    title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
+                    Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     self.xlabel = "Range (Km)"
                     self.ylabel = "Intensity - IQ"
                     self.y = yreal
                     self.x = x
-                    self.titles[0] = title
                     for i,ax in enumerate(self.axes):
                         title = "Channel %d" %(i)
                         if ax.firsttime:
                         else:
                             ax.plt_r.set_data(x, yreal[i,:])
                             ax.plt_i.set_data(x, yimag[i,:])
+                    plt.suptitle(Maintitle)
                 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
                     y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
                         else:
                             ax.plt_r.set_data(x, ychannel)
                 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):

schainpy/model/io/jroIO_base.py +84 -84

             def isRadarFile(file):
                 try:
                     year = int(file[1:5])
                     doy = int(file[5:8])
                     set = int(file[8:11])
             def getDateFromRadarFile(file):
                 try:
                     year = int(file[1:5])
                     doy = int(file[5:8])
                     set = int(file[8:11])
                 except:
                     return None
                 return thisDate
             def parse_format(s, fmt):
                 for i in range(fmt.count('%')):
                     x = fmt.index('%')
                     d = DT_DIRECTIVES[fmt[x:x+2]]
                 def run(self):
                     raise NotImplementedError
                 def getAllowedArgs(self):
                     if hasattr(self, '__attrs__'):
                     for key, value in kwargs.items():
                         setattr(self, key, value)
                 def find_folders(self, path, startDate, endDate, folderfmt, last=False):
                     folders = [x for f in path.split(',')
                                for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
                     folders.sort()
                     if last:
                         folders = [folders[-1]]
                     for folder in folders:
                         try:
                             dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
                             if dt >= startDate and dt <= endDate:
                                 yield os.path.join(path, folder)
                             else:
                             log.log('Skiping folder {}'.format(folder), self.name)
                             continue
                     return
                 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
                                expLabel='', last=False):
                     for path in folders:
                         files = glob.glob1(path, '*{}'.format(ext))
                         files.sort()
                         if last:
                             if files:
                                 fo = files[-1]
                                 try:
                                     dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
                                     yield os.path.join(path, expLabel, fo)
                                 except Exception as e:
                                     pass
                                 return
                             else:
                                 return
                         for fo in files:
                             try:
                                 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
                                 if dt >= startDate and dt <= endDate:
                                     yield os.path.join(path, expLabel, fo)
                                 else:
                                     log.log('Skiping file {}'.format(fo), self.name)
                             except Exception as e:
                                 log.log('Skiping file {}'.format(fo), self.name)
                                 continue
                 def searchFilesOffLine(self, path, startDate, endDate,
                                        expLabel, ext, walk,
                                        filefmt, folderfmt):
                     """Search files in offline mode for the given arguments
                             path, startDate, endDate, folderfmt)
                     else:
                         folders = path.split(',')
                     return self.find_files(
                         folders, ext, filefmt, startDate, endDate, expLabel)
                 def searchFilesOnLine(self, path, startDate, endDate,
                                       expLabel, ext, walk,
                                       filefmt, folderfmt):
                     """Search for the last file of the last folder
                     Return:
                         generator with the full path of last filename
                     """
                     if walk:
                         folders = self.find_folders(
                             path, startDate, endDate, folderfmt, last=True)
                     else:
                         folders = path.split(',')
                     return self.find_files(
                         folders, ext, filefmt, startDate, endDate, expLabel, last=True)
                     while True:
                         if self.fp != None:
                             self.fp.close()
                         if self.online:
                             newFile = self.setNextFileOnline()
                         else:
                             newFile = self.setNextFileOffline()
                         if not(newFile):
                             if self.online:
                                 raise schainpy.admin.SchainError('Time to wait for new files reach')
                                     raise schainpy.admin.SchainWarning('No files found in the given path')
                                 else:
                                     raise schainpy.admin.SchainWarning('No more files to read')
                         if self.verifyFile(self.filename):
                             break
                     log.log('Opening file: %s' % self.filename, self.name)
                     self.readFirstHeader()
                         self.filename
                         self.fp
                         self.filesize
                     Return:
                         boolean
                     nextFile = True
                     nextDay = False
                     for nFiles in range(self.nFiles+1):
                         for nTries in range(self.nTries):
                             fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
                             if fullfilename is not None:
                                 self.name)
                             time.sleep(self.delay)
                             nextFile = False
                             continue
                         if fullfilename is not None:
                             break
                         self.nTries = 1
                         nextFile = True
                         if nFiles == (self.nFiles - 1):
                             log.log('Trying with next day...', self.name)
                             nextDay = True
                             self.nTries = 3
                     if fullfilename:
                         self.fileSize = os.path.getsize(fullfilename)
                         self.flagNoMoreFiles = 0
                         self.fileIndex += 1
                         return 1
                     else:
                         return 0
                 def setNextFileOffline(self):
                     """Open the next file to be readed in offline mode"""
                     try:
                         filename = next(self.filenameList)
                         self.fileIndex +=1
                     except StopIteration:
                         self.flagNoMoreFiles = 1
                         return 0
                     self.filename = filename
                     self.fileSize = os.path.getsize(filename)
                     self.flagIsNewFile = 1
                     return 1
                 @staticmethod
                 def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
                     """Check if the given datetime is in range"""
+                    startDateTime= datetime.datetime.combine(startDate,startTime)
-                    if startDate <= dt.date() <= endDate:
+                    endDateTime = datetime.datetime.combine(endDate,endTime)
-                        if startTime <= dt.time() <= endTime:
+                    if startDateTime <= dt <= endDateTime:
-                            return True
+                        return True
                     return False
                 def verifyFile(self, filename):
                     """Check for a valid file
                     Arguments:
                         filename -- full path filename
                     Return:
                         boolean
                     """
                     """Check if the next file to be readed exists"""
                     raise NotImplementedError
                 def readFirstHeader(self):
                     """Parse the file header"""
                     Return:
                         str -- fullpath of the file
                     """
                     if nextFile:
                         self.set += 1
                     if nextDay:
                         prefixFileList = ['d', 'D']
                     elif self.ext.lower() == ".pdata":  # spectra
                         prefixFileList = ['p', 'P']
                     # barrido por las combinaciones posibles
                     for prefixDir in prefixDirList:
                         thispath = self.path
                             if os.path.exists(fullfilename):
                                 return fullfilename, filename
                     return None, filename
                 def __waitNewBlock(self):
                     """
                     Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
                 def __setNewBlock(self):
                     if self.fp == None:
                         return 0
                     if self.flagIsNewFile:
                         self.lastUTTime = self.basicHeaderObj.utc
                         return 1
                     currentSize = self.fileSize - self.fp.tell()
                     neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
                     if (currentSize >= neededSize):
                         self.basicHeaderObj.read(self.fp)
                         self.lastUTTime = self.basicHeaderObj.utc
                         return 1
                     if self.__waitNewBlock():
                         self.lastUTTime = self.basicHeaderObj.utc
                         return 1
                     except IOError:
                         log.error("File {} can't be opened".format(filename), self.name)
                         return False
                     if self.online and self.waitDataBlock(0):
                         pass
                     basicHeaderObj = BasicHeader(LOCALTIME)
                     systemHeaderObj = SystemHeader()
                     radarControllerHeaderObj = RadarControllerHeader()
                         dt2 = basicHeaderObj.datatime
                         if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
                                 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
                             flag = False
                     fp.close()
                     return flag
                     return dateList
                 def setup(self, **kwargs):
                     self.set_kwargs(**kwargs)
                     if not self.ext.startswith('.'):
                         self.ext = '.{}'.format(self.ext)
                     if self.server is not None:
                         if 'tcp://' in self.server:
                             address = server
                             for nTries in range(self.nTries):
                                 fullpath = self.searchFilesOnLine(self.path, self.startDate,
                                     self.endDate, self.expLabel, self.ext, self.walk,
                                     self.filefmt, self.folderfmt)
                                 try:
                                     fullpath = next(fullpath)
                                 except:
                                     fullpath = None
                                 if fullpath:
                                     break
                                 log.warning(
                                     'Waiting {} sec for a valid file in {}: try {} ...'.format(
                                         self.delay, self.path, nTries + 1),
                                     self.name)
                                 time.sleep(self.delay)
                             if not(fullpath):
                                 raise schainpy.admin.SchainError(
                                     'There isn\'t any valid file in {}'.format(self.path))
                             pathname, filename = os.path.split(fullpath)
                             self.year = int(filename[1:5])
                             self.doy = int(filename[5:8])
                             self.set = int(filename[8:11]) - 1
                         else:
                             log.log("Searching files in {}".format(self.path), self.name)
                             self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
                                 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
                         self.setNextFile()
                     return
                     self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
                     self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
                 def getFirstHeader(self):
                     raise NotImplementedError
                     """
                     Arguments:
                         path        :
                         startDate   :
                         endDate     :
                         startTime   :
                         endTime     :
                     dtype_width = get_dtype_width(dtype_index)
                     return dtype_width
                 def getProcessFlags(self):
                     processFlags = 0
                     self.basicHeaderObj.size = self.basicHeaderSize  # bytes
                     self.basicHeaderObj.version = self.versionFile
                     self.basicHeaderObj.dataBlock = self.nTotalBlocks
                     utc = numpy.floor(self.dataOut.utctime)
                     milisecond = (self.dataOut.utctime - utc) * 1000.0
                     self.basicHeaderObj.utc = utc
                     self.basicHeaderObj.miliSecond = milisecond
                     self.basicHeaderObj.timeZone = self.dataOut.timeZone
                     if self.dataOut.datatime.date() > self.fileDate:
                         setFile = 0
                         self.nTotalBlocks = 0
                     filen = '{}{:04d}{:03d}{:03d}{}'.format(
                             self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
                     filename = os.path.join(path, subfolder, filen)
                     self.ext = ext.lower()
                     self.path = path
                     if set is None:
                         self.setFile = -1
                     else:
                         self.setFile = set - 1
                     self.blocksPerFile = blocksPerFile
                     self.profilesPerBlock = profilesPerBlock

schainpy/model/io/jroIO_kamisr.py +5 -3

                     else:
                         #get the last file - 1
-                        self.filenameList = [self.filenameList[-1]]
+                        self.filenameList = [self.filenameList[-2]]
                     new_dirnameList = []
                     for dirname in self.dirnameList:
                         junk = numpy.array([dirname in x for x in self.filenameList])
                         self.__selectDataForTimes(online=True)
                         filename = self.filenameList[0]
                         wait = 0
-                        #self.__waitForNewFile=180 ## DEBUG:
+                        self.__waitForNewFile=300 ## DEBUG:
                         while self.__filename_online == filename:
                             print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
                             if wait == 5:
                     #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
                     #print(self.dataOut.utctime)
                     self.dataOut.profileIndex = self.profileIndex
+                    #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
                     self.dataOut.flagNoData = False
             #         if indexprof == 0:
             #             print self.dataOut.utctime
                     self.profileIndex += 1
-                    return self.dataOut.data
+                    #return self.dataOut.data
                 def run(self, **kwargs):
                         self.isConfig = True
                     self.getData()
+                    return(self.dataOut.data)

schainpy/model/proc/jroproc_heispectra.py +3 -4

             from schainpy.utils import log
             class SpectraHeisProc(ProcessingUnit):
                 def __init__(self):#, **kwargs):
                     if self.dataIn.type == "Fits":
                         self.__updateObjFromFits()
                         self.dataOut.flagNoData = False
                         return
                     if self.dataIn.type == "SpectraHeis":
                         self.dataOut.copy(self.dataIn)
             #             dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
             #            dataOut.timeInterval = self.__timeInterval*self.n
                         dataOut.flagNoData = False
-                    return dataOut
  No newline at end of file
+                    return dataOut

schainpy/model/proc/jroproc_spectra.py +559 0

This diff has been collapsed as it changes many lines, (559 lines changed) Show them Hide them
	@@ -12,12 +12,15 import time
	12	import itertools	12	import itertools
	13		13
	14	import numpy	14	import numpy
			15	import math
	15		16
	16	from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation	17	from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
	17	from schainpy.model.data.jrodata import Spectra	18	from schainpy.model.data.jrodata import Spectra
	18	from schainpy.model.data.jrodata import hildebrand_sekhon	19	from schainpy.model.data.jrodata import hildebrand_sekhon
	19	from schainpy.utils import log	20	from schainpy.utils import log
	20		21
			22	from scipy.optimize import curve_fit
			23
	21		24
	22	class SpectraProc(ProcessingUnit):	25	class SpectraProc(ProcessingUnit):
	23		26
	@@ -478,6 +481,562 class removeDC(Operation):
	478		481
	479	return self.dataOut	482	return self.dataOut
	480		483
			484	# import matplotlib.pyplot as plt
			485
			486	def fit_func( x, a0, a1, a2): #, a3, a4, a5):
			487	z = (x - a1) / a2
			488	y = a0 * numpy.exp(-z*2 / a2) #+ a3 + a4 x + a5 * x**2
			489	return y
			490	class CleanRayleigh(Operation):
			491
			492	def __init__(self):
			493
			494	Operation.__init__(self)
			495	self.i=0
			496	self.isConfig = False
			497	self.__dataReady = False
			498	self.__profIndex = 0
			499	self.byTime = False
			500	self.byProfiles = False
			501
			502	self.bloques = None
			503	self.bloque0 = None
			504
			505	self.index = 0
			506
			507	self.buffer = 0
			508	self.buffer2 = 0
			509	self.buffer3 = 0
			510	#self.min_hei = None
			511	#self.max_hei = None
			512
			513	def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
			514
			515	self.nChannels = dataOut.nChannels
			516	self.nProf = dataOut.nProfiles
			517	self.nPairs = dataOut.data_cspc.shape[0]
			518	self.pairsArray = numpy.array(dataOut.pairsList)
			519	self.spectra = dataOut.data_spc
			520	self.cspectra = dataOut.data_cspc
			521	self.heights = dataOut.heightList #alturas totales
			522	self.nHeights = len(self.heights)
			523	self.min_hei = min_hei
			524	self.max_hei = max_hei
			525	if (self.min_hei == None):
			526	self.min_hei = 0
			527	if (self.max_hei == None):
			528	self.max_hei = dataOut.heightList[-1]
			529	self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
			530	self.heightsClean = self.heights[self.hval] #alturas filtradas
			531	self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
			532	self.nHeightsClean = len(self.heightsClean)
			533	self.channels = dataOut.channelList
			534	self.nChan = len(self.channels)
			535	self.nIncohInt = dataOut.nIncohInt
			536	self.__initime = dataOut.utctime
			537	self.maxAltInd = self.hval[-1]+1
			538	self.minAltInd = self.hval[0]
			539
			540	self.crosspairs = dataOut.pairsList
			541	self.nPairs = len(self.crosspairs)
			542	self.normFactor = dataOut.normFactor
			543	self.nFFTPoints = dataOut.nFFTPoints
			544	self.ippSeconds = dataOut.ippSeconds
			545	self.currentTime = self.__initime
			546	self.pairsArray = numpy.array(dataOut.pairsList)
			547	self.factor_stdv = factor_stdv
			548
			549
			550
			551	if n != None :
			552	self.byProfiles = True
			553	self.nIntProfiles = n
			554	else:
			555	self.__integrationtime = timeInterval
			556
			557	self.__dataReady = False
			558	self.isConfig = True
			559
			560
			561
			562	def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
			563	#print (dataOut.utctime)
			564	if not self.isConfig :
			565	#print("Setting config")
			566	self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
			567	#print("Config Done")
			568	tini=dataOut.utctime
			569
			570	if self.byProfiles:
			571	if self.__profIndex == self.nIntProfiles:
			572	self.__dataReady = True
			573	else:
			574	if (tini - self.__initime) >= self.__integrationtime:
			575	#print(tini - self.__initime,self.__profIndex)
			576	self.__dataReady = True
			577	self.__initime = tini
			578
			579	#if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
			580
			581	if self.__dataReady:
			582	print("Data ready",self.__profIndex)
			583	self.__profIndex = 0
			584	jspc = self.buffer
			585	jcspc = self.buffer2
			586	#jnoise = self.buffer3
			587	self.buffer = dataOut.data_spc
			588	self.buffer2 = dataOut.data_cspc
			589	#self.buffer3 = dataOut.noise
			590	self.currentTime = dataOut.utctime
			591	if numpy.any(jspc) :
			592	#print( jspc.shape, jcspc.shape)
			593	jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
			594	jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
			595	self.__dataReady = False
			596	#print( jspc.shape, jcspc.shape)
			597	dataOut.flagNoData = False
			598	else:
			599	dataOut.flagNoData = True
			600	self.__dataReady = False
			601	return dataOut
			602	else:
			603	#print( len(self.buffer))
			604	if numpy.any(self.buffer):
			605	self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
			606	self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
			607	self.buffer3 += dataOut.data_dc
			608	else:
			609	self.buffer = dataOut.data_spc
			610	self.buffer2 = dataOut.data_cspc
			611	self.buffer3 = dataOut.data_dc
			612	#print self.index, self.fint
			613	#print self.buffer2.shape
			614	dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
			615	self.__profIndex += 1
			616	return dataOut ## NOTE: REV
			617
			618	# if self.index == 0 and self.fint == 1 :
			619	# if jspc != None:
			620	# print len(jspc), jspc.shape
			621	# jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4))
			622	# print jspc.shape
			623	# dataOut.flagNoData = True
			624	# return dataOut
			625	# if path != None:
			626	# sys.path.append(path)
			627	# self.library = importlib.import_module(file)
			628	#
			629	# To be inserted as a parameter
			630
			631	#Set constants
			632	#constants = self.library.setConstants(dataOut)
			633	#dataOut.constants = constants
			634
			635	#snrth= 20
			636
			637	#crosspairs = dataOut.groupList
			638	#noise = dataOut.noise
			639	#print( nProf,heights)
			640	#print( jspc.shape, jspc.shape[0])
			641	#print noise
			642	#print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise)
			643	#jnoise = jnoise/N
			644	#noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
			645	#print( noise)
			646	#power = numpy.sum(spectra, axis=1)
			647	#print power[0,:]
			648	#print("CROSSPAIRS",crosspairs)
			649	#nPairs = len(crosspairs)
			650	#print(numpy.shape(dataOut.data_spc))
			651
			652	#absc = dataOut.abscissaList[:-1]
			653
			654	#print absc.shape
			655	#nBlocks=149
			656	#print('spectra', spectra.shape)
			657	#print('noise print', crosspairs)
			658	#print('spectra', spectra.shape)
			659	#print('cspectra', cspectra.shape)
			660	#print numpy.array(dataOut.data_pre[1]).shape
			661	#spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs)
			662
			663
			664	#index = tini.tm_hour*12+tini.tm_min/5
			665	# jspc = jspc/self.nFFTPoints/self.normFactor
			666	# jcspc = jcspc/self.nFFTPoints/self.normFactor
			667
			668
			669
			670	#dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
			671	#tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei)
			672	tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
			673	#jspectra = tmp_spectra*len(jspc[:,0,0,0])
			674	#jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
			675
			676	dataOut.data_spc = tmp_spectra
			677	dataOut.data_cspc = tmp_cspectra
			678	dataOut.data_dc = self.buffer3
			679	dataOut.nIncohInt *= self.nIntProfiles
			680	dataOut.utctime = self.currentTime #tiempo promediado
			681	print("Time: ",time.localtime(dataOut.utctime))
			682	# dataOut.data_spc = sat_spectra
			683	# dataOut.data_cspc = sat_cspectra
			684	self.buffer = 0
			685	self.buffer2 = 0
			686	self.buffer3 = 0
			687
			688	return dataOut
			689
			690	def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
			691	print("OP cleanRayleigh")
			692	#import matplotlib.pyplot as plt
			693	#for k in range(149):
			694
			695	rfunc = cspectra.copy() #self.bloques
			696	val_spc = spectra0.0 #self.bloque00.0
			697	val_cspc = cspectra0.0 #self.bloques0.0
			698	in_sat_spectra = spectra.copy() #self.bloque0
			699	in_sat_cspectra = cspectra.copy() #self.bloques
			700
			701	raxs = math.ceil(math.sqrt(self.nPairs))
			702	caxs = math.ceil(self.nPairs/raxs)
			703
			704	#print(self.hval)
			705	#print numpy.absolute(rfunc[:,0,0,14])
			706	for ih in range(self.minAltInd,self.maxAltInd):
			707	for ifreq in range(self.nFFTPoints):
			708	# fig, axs = plt.subplots(raxs, caxs)
			709	# fig2, axs2 = plt.subplots(raxs, caxs)
			710	col_ax = 0
			711	row_ax = 0
			712	for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
			713	#print("ii: ",ii)
			714	if (col_ax%caxs==0 and col_ax!=0):
			715	col_ax = 0
			716	row_ax += 1
			717	func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
			718	#print(func2clean.shape)
			719	val = (numpy.isfinite(func2clean)==True).nonzero()
			720
			721	if len(val)>0: #limitador
			722	min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
			723	if min_val <= -40 :
			724	min_val = -40
			725	max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
			726	if max_val >= 200 :
			727	max_val = 200
			728	#print min_val, max_val
			729	step = 1
			730	#print("Getting bins and the histogram")
			731	x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
			732	y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
			733	#print(len(y_dist),len(binstep[:-1]))
			734	#print(row_ax,col_ax, " ..")
			735	#print(self.pairsArray[ii][0],self.pairsArray[ii][1])
			736	mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
			737	sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
			738	parg = [numpy.amax(y_dist),mean,sigma]
			739	gauss_fit, covariance = None, None
			740	newY = None
			741	try :
			742	gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
			743	mode = gauss_fit[1]
			744	stdv = gauss_fit[2]
			745	#print(" FIT OK",gauss_fit)
			746	'''
			747	newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
			748	axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
			749	axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
			750	axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
			751	except:
			752	mode = mean
			753	stdv = sigma
			754	#print("FIT FAIL")
			755
			756
			757	#print(mode,stdv)
			758	#Removing echoes greater than mode + 3*stdv
			759	#factor_stdv = 2
			760	noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
			761	#noval tiene los indices que se van a remover
			762	#print("Pair ",ii," novals: ",len(noval[0]))
			763	if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
			764	novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
			765	#print(novall)
			766	#print(" ")
			767	cross_pairs = self.pairsArray[ii]
			768	#Getting coherent echoes which are removed.
			769	if len(novall[0]) > 0:
			770
			771	val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
			772	val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
			773	val_cspc[novall[0],ii,ifreq,ih] = 1
			774	#print("OUT NOVALL 1")
			775	#Removing coherent from ISR data
			776
			777	#print(spectra[:,ii,ifreq,ih])
			778	new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
			779	mean_cspc = numpy.mean(new_a)
			780	new_b = numpy.delete(spectra[:,cross_pairs[0],ifreq,ih], noval[0])
			781	mean_spc0 = numpy.mean(new_b)
			782	new_c = numpy.delete(spectra[:,cross_pairs[1],ifreq,ih], noval[0])
			783	mean_spc1 = numpy.mean(new_c)
			784	spectra[noval,cross_pairs[0],ifreq,ih] = mean_spc0
			785	spectra[noval,cross_pairs[1],ifreq,ih] = mean_spc1
			786	cspectra[noval,ii,ifreq,ih] = mean_cspc
			787
			788	'''
			789	func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
			790	y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
			791	axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
			792	axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
			793	axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
			794	'''
			795
			796	col_ax += 1 #contador de ploteo columnas
			797	##print(col_ax)
			798	'''
			799	title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
			800	title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
			801	fig.suptitle(title)
			802	fig2.suptitle(title2)
			803	plt.show()'''
			804
			805	''' channels = channels
			806	cross_pairs = cross_pairs
			807	#print("OUT NOVALL 2")
			808
			809	vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
			810	vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
			811	vcross = numpy.concatenate((vcross0,vcross1),axis=None)
			812	#print('vcros =', vcross)
			813
			814	#Getting coherent echoes which are removed.
			815	if len(novall) > 0:
			816	#val_spc[novall,ii,ifreq,ih] = 1
			817	val_spc[ii,ifreq,ih,novall] = 1
			818	if len(vcross) > 0:
			819	val_cspc[vcross,ifreq,ih,novall] = 1
			820
			821	#Removing coherent from ISR data.
			822	self.bloque0[ii,ifreq,ih,noval] = numpy.nan
			823	if len(vcross) > 0:
			824	self.bloques[vcross,ifreq,ih,noval] = numpy.nan
			825	'''
			826
			827	print("Getting average of the spectra and cross-spectra from incoherent echoes.")
			828	out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
			829	out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
			830	for ih in range(self.nHeights):
			831	for ifreq in range(self.nFFTPoints):
			832	for ich in range(self.nChan):
			833	tmp = spectra[:,ich,ifreq,ih]
			834	valid = (numpy.isfinite(tmp[:])==True).nonzero()
			835	# if ich == 0 and ifreq == 0 and ih == 17 :
			836	# print tmp
			837	# print valid
			838	# print len(valid[0])
			839	#print('TMP',tmp)
			840	if len(valid[0]) >0 :
			841	out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
			842	#for icr in range(nPairs):
			843	for icr in range(self.nPairs):
			844	tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
			845	valid = (numpy.isfinite(tmp)==True).nonzero()
			846	if len(valid[0]) > 0:
			847	out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
			848	'''
			849	# print('##########################################################')
			850	print("Removing fake coherent echoes (at least 4 points around the point)")
			851
			852	val_spectra = numpy.sum(val_spc,0)
			853	val_cspectra = numpy.sum(val_cspc,0)
			854
			855	val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
			856	val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
			857
			858	for i in range(nChan):
			859	for j in range(nProf):
			860	for k in range(nHeights):
			861	if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
			862	val_spc[:,i,j,k] = 0.0
			863	for i in range(nPairs):
			864	for j in range(nProf):
			865	for k in range(nHeights):
			866	if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
			867	val_cspc[:,i,j,k] = 0.0
			868
			869	# val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProfnHeightsnChan))
			870	# if numpy.isfinite(val_spectra)==str(True):
			871	# noval = (val_spectra<1).nonzero()
			872	# if len(noval) > 0:
			873	# val_spc[:,noval] = 0.0
			874	# val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights))
			875
			876	#val_cspc = numpy.reshape(val_spc, (149,nChannHeightsnProf))
			877	#if numpy.isfinite(val_cspectra)==str(True):
			878	# noval = (val_cspectra<1).nonzero()
			879	# if len(noval) > 0:
			880	# val_cspc[:,noval] = 0.0
			881	# val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights))
			882	tmp_sat_spectra = spectra.copy()
			883	tmp_sat_spectra = tmp_sat_spectra*numpy.nan
			884	tmp_sat_cspectra = cspectra.copy()
			885	tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
			886	'''
			887	# fig = plt.figure(figsize=(6,5))
			888	# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
			889	# ax = fig.add_axes([left, bottom, width, height])
			890	# cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
			891	# ax.clabel(cp, inline=True,fontsize=10)
			892	# plt.show()
			893	'''
			894	val = (val_spc > 0).nonzero()
			895	if len(val[0]) > 0:
			896	tmp_sat_spectra[val] = in_sat_spectra[val]
			897	val = (val_cspc > 0).nonzero()
			898	if len(val[0]) > 0:
			899	tmp_sat_cspectra[val] = in_sat_cspectra[val]
			900
			901	print("Getting average of the spectra and cross-spectra from incoherent echoes 2")
			902	sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float)
			903	sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex)
			904	for ih in range(nHeights):
			905	for ifreq in range(nProf):
			906	for ich in range(nChan):
			907	tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
			908	valid = (numpy.isfinite(tmp)).nonzero()
			909	if len(valid[0]) > 0:
			910	sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
			911
			912	for icr in range(nPairs):
			913	tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
			914	valid = (numpy.isfinite(tmp)).nonzero()
			915	if len(valid[0]) > 0:
			916	sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
			917	'''
			918	#self.__dataReady= True
			919	#sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
			920	#if not self.__dataReady:
			921	#return None, None
			922	#return out_spectra, out_cspectra ,sat_spectra,sat_cspectra
			923	return out_spectra, out_cspectra
			924
			925	def REM_ISOLATED_POINTS(self,array,rth):
			926	# import matplotlib.pyplot as plt
			927	if rth == None :
			928	rth = 4
			929	print("REM ISO")
			930	num_prof = len(array[0,:,0])
			931	num_hei = len(array[0,0,:])
			932	n2d = len(array[:,0,0])
			933
			934	for ii in range(n2d) :
			935	#print ii,n2d
			936	tmp = array[ii,:,:]
			937	#print tmp.shape, array[ii,101,:],array[ii,102,:]
			938
			939	# fig = plt.figure(figsize=(6,5))
			940	# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
			941	# ax = fig.add_axes([left, bottom, width, height])
			942	# x = range(num_prof)
			943	# y = range(num_hei)
			944	# cp = ax.contour(y,x,tmp)
			945	# ax.clabel(cp, inline=True,fontsize=10)
			946	# plt.show()
			947
			948	#indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
			949	tmp = numpy.reshape(tmp,num_prof*num_hei)
			950	indxs1 = (numpy.isfinite(tmp)==True).nonzero()
			951	indxs2 = (tmp > 0).nonzero()
			952
			953	indxs1 = (indxs1[0])
			954	indxs2 = indxs2[0]
			955	#indxs1 = numpy.array(indxs1[0])
			956	#indxs2 = numpy.array(indxs2[0])
			957	indxs = None
			958	#print indxs1 , indxs2
			959	for iv in range(len(indxs2)):
			960	indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
			961	#print len(indxs2), indv
			962	if len(indv[0]) > 0 :
			963	indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
			964	# print indxs
			965	indxs = indxs[1:]
			966	#print(indxs, len(indxs))
			967	if len(indxs) < 4 :
			968	array[ii,:,:] = 0.
			969	return
			970
			971	xpos = numpy.mod(indxs ,num_hei)
			972	ypos = (indxs / num_hei)
			973	sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
			974	#print sx
			975	xpos = xpos[sx]
			976	ypos = ypos[sx]
			977
			978	# ********************************* Cleaning isolated points ********************************
			979	ic = 0
			980	while True :
			981	r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
			982	#no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
			983	#plt.plot(r)
			984	#plt.show()
			985	no_coh1 = (numpy.isfinite(r)==True).nonzero()
			986	no_coh2 = (r <= rth).nonzero()
			987	#print r, no_coh1, no_coh2
			988	no_coh1 = numpy.array(no_coh1[0])
			989	no_coh2 = numpy.array(no_coh2[0])
			990	no_coh = None
			991	#print valid1 , valid2
			992	for iv in range(len(no_coh2)):
			993	indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
			994	if len(indv[0]) > 0 :
			995	no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
			996	no_coh = no_coh[1:]
			997	#print len(no_coh), no_coh
			998	if len(no_coh) < 4 :
			999	#print xpos[ic], ypos[ic], ic
			1000	# plt.plot(r)
			1001	# plt.show()
			1002	xpos[ic] = numpy.nan
			1003	ypos[ic] = numpy.nan
			1004
			1005	ic = ic + 1
			1006	if (ic == len(indxs)) :
			1007	break
			1008	#print( xpos, ypos)
			1009
			1010	indxs = (numpy.isfinite(list(xpos))==True).nonzero()
			1011	#print indxs[0]
			1012	if len(indxs[0]) < 4 :
			1013	array[ii,:,:] = 0.
			1014	return
			1015
			1016	xpos = xpos[indxs[0]]
			1017	ypos = ypos[indxs[0]]
			1018	for i in range(0,len(ypos)):
			1019	ypos[i]=int(ypos[i])
			1020	junk = tmp
			1021	tmp = junk*0.0
			1022
			1023	tmp[list(xpos + (yposnum_hei))] = junk[list(xpos + (yposnum_hei))]
			1024	array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
			1025
			1026	#print array.shape
			1027	#tmp = numpy.reshape(tmp,(num_prof,num_hei))
			1028	#print tmp.shape
			1029
			1030	# fig = plt.figure(figsize=(6,5))
			1031	# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
			1032	# ax = fig.add_axes([left, bottom, width, height])
			1033	# x = range(num_prof)
			1034	# y = range(num_hei)
			1035	# cp = ax.contour(y,x,array[ii,:,:])
			1036	# ax.clabel(cp, inline=True,fontsize=10)
			1037	# plt.show()
			1038	return array
			1039
	481	class removeInterference(Operation):	1040	class removeInterference(Operation):
	482		1041
	483	def removeInterference2(self):	1042	def removeInterference2(self):

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages