schain Commit - r1548:334daf8eec15 · Jicamarca Repository

remSatsProfiles

joabAM -

r1548:334daf8eec15

parent child

Context file:

r1548:334daf8eec15

Collapse all files

schainpy/model/data/jrodata.py +10 -5

@@ -606,11 +606,16 class Spectra(JROData):
606	def getPower(self):	606	def getPower(self):
607		607
608	factor = self.normFactor	608	factor = self.normFactor
609	z = numpy.divide(self.data_spc,factor)	609	power = numpy.zeros( (self.nChannels,self.nHeights) )
610	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)	610	for ch in range(self.nChannels):
611	avg = numpy.average(z, axis=1)	611	if hasattr(factor,'shape'):
612		612	z = numpy.divide(self.data_spc[ch],factor[ch])
613	return 10 * numpy.log10(avg)	613	else:
		614	z = numpy.divide(self.data_spc[ch],factor)
		615	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
		616	avg = numpy.average(z, axis=0)
		617	power[ch,:] = 10 * numpy.log10(avg)
		618	return power
614		619
615	def getCoherence(self, pairsList=None, phase=False):	620	def getCoherence(self, pairsList=None, phase=False):
616		621

schainpy/model/graphics/jroplot_spectra.py +1 -1

                     data = {}
                     meta = {}
                     norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt  * dataOut.windowOfFilter
-                    #noise =  10*numpy.log10(dataOut.getNoise()/norm)
                     noise =  10*numpy.log10(dataOut.getNoise())
                     noise = noise.reshape(dataOut.nChannels, 1)
                     data['noise'] = noise
                     data['noise'] = n0
                     noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
                     data['noiseless_rti'] = dataOut.getPower() - noise
                     return data, meta

schainpy/model/proc/jroproc_parameters.py +33 -12

                     #----------------------    Spectra Data    ---------------------------
                     if self.dataIn.type == "Spectra":
                         self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
                         self.dataOut.data_spc = self.dataIn.data_spc
                         self.dataOut.data_cspc = self.dataIn.data_cspc
                             self.dataOut.groupList = self.dataIn.pairsList
                         self.dataOut.flagNoData = False
+                        self.dataOut.noise_estimation = self.dataIn.noise_estimation
                         if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
                             self.dataOut.ChanDist = self.dataIn.ChanDist
                         else: self.dataOut.ChanDist = None
                     nChannel = data.shape[0]
                     data_param = numpy.zeros((nChannel, 4, data.shape[2]))
+                    #norm = dataOut.nIncohInt/dataOut.max_nIncohInt
+                    norm = dataOut.max_nIncohInt/dataOut.nIncohInt
                     for ind in range(nChannel):
-                        data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
+                        if hasattr(norm, "shape"):
+                            data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind], normFactor=norm[ind] )
+                        else:
+                            data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
                     dataOut.moments = data_param[:,1:,:]
                     dataOut.data_snr = data_param[:,0]
                     dataOut.data_dop = data_param[:,2]
                     dataOut.data_width = data_param[:,3]
                     return dataOut
-                def __calculateMoments(self, oldspec, oldfreq, n0,
+                def __calculateMoments(self, oldspec, oldfreq, n0,normFactor = 1,
                                        nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
                     if (nicoh is None): nicoh = 1
                     if (n0 < 1.e-20):   n0 = 1.e-20
                     freq = oldfreq
-                    vec_power = numpy.zeros(oldspec.shape[1])
+                    nheights = oldspec.shape[1]
-                    vec_fd = numpy.zeros(oldspec.shape[1])
+                    vec_power = numpy.zeros(nheights)
-                    vec_w = numpy.zeros(oldspec.shape[1])
+                    vec_fd = numpy.zeros(nheights)
-                    vec_snr = numpy.zeros(oldspec.shape[1])
+                    vec_w = numpy.zeros(nheights)
+                    vec_snr = numpy.zeros(nheights)
+                    norm = 1
                     # oldspec = numpy.ma.masked_invalid(oldspec)
-                    for ind in range(oldspec.shape[1]):
+                    for ind in range(nheights):
                         spec = oldspec[:,ind]
                         aux = spec*fwindow
                         max_spec = aux.max()
                         m = aux.tolist().index(max_spec)
+                        if hasattr(normFactor, "shape"):
+                            norm = normFactor[ind]
                         # Smooth
                         if (smooth == 0):
                             spec2 = spec
                             spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
                         # Moments Estimation
                         bb = spec2[numpy.arange(m,spec2.size)]
                         bb = (bb<n0).nonzero()
                         bb = bb[0]
                         if (ss1 > m):
                             ss1 = m
                         valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
                         signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean()    # D. Scipión added with correct definition
                         power = ((spec2[valid] - n0) * fwindow[valid]).sum()
                         fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
                         w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
-                        snr = (spec2.mean()-n0)/n0
-                        if (snr < 1.e-20) :
+                        spec2 *= (norm*norm)   #compensation for sats remove
-                            snr = 1.e-20
+                        #snr = (spec2.mean()-n0)/n0
+                        snr = (spec2.mean())/n0
+                        #print(snr)
+                        # if  snr < 1.0:  # By definition SNR must be higher than 1
+                        #     snr = 0
                         # vec_power[ind] = power    #D. Scipión replaced with the line below
                         vec_power[ind] = total_power

schainpy/model/proc/jroproc_spectra.py +12 -11

                     else:
                         raise ValueError("The type of input object {} is not valid".format(
                             self.dataIn.type))
+                    #print("spc proc Done")
                 def __selectPairs(self, pairsList):
                         noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
                     self.dataOut.noise_estimation = noise.copy() # dataOut.noise
                     #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
+                    #print("2: ",self.dataOut.noise_estimation)
                     #print(self.dataOut.flagNoData)
                     #print("getNoise Done", noise)
                     return self.dataOut
                     #self.ByLags = dataOut.ByLags   ###REDEFINIR
                     self.ByLags = False
                     self.maxProfilesInt = 1
+                    self.__nChannels = dataOut.nChannels
                     if DPL != None:
                         self.DPL=DPL
                     else:
                     ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
                     self.minHei_ind = ind_list1[0][0]
                     self.maxHei_ind = ind_list2[0][-1]
+                    #print("setup rem sats done")
                 def putData(self, data_spc, data_cspc, data_dc):
                     """
                     self.__buffer_spc.append(data_spc)
-                    if data_cspc is None:
+                    if self.__nChannels < 2:
                         self.__buffer_cspc = None
                     else:
                         self.__buffer_cspc.append(data_cspc)
                     buffer_cspc=None
                     #print("aes: ", self.__buffer_cspc)
                     self.__buffer_spc=numpy.array(self.__buffer_spc)
-                    if self.__buffer_cspc != None:
+                    if self.__nChannels > 1 :
                         self.__buffer_cspc=numpy.array(self.__buffer_cspc)
                     #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
                     self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
                     for k in range(self.minHei_ind,self.maxHei_ind):
-                        if self.__buffer_cspc != None:
+                        if self.__nChannels > 1:
                             buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
                         outliers_IDs_cspc=[]
                             # plt.show()
                             if indexmin != buffer1.shape[0]:
-                                if self.nChannels > 1:
+                                if self.__nChannels > 1:
                                     cspc_outliers_exist= True
                                 lt=outliers_IDs
                             self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
-                            if self.__buffer_cspc != None:
+                            if self.__nChannels > 1:
                                 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
-                        if self.__buffer_cspc != None:
+                        if self.__nChannels > 1:
                             outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
                         if cspc_outliers_exist:
                                 #buffer_cspc[p,:]=avg
                                 buffer_cspc[p,:] = numpy.NaN
-                        if self.__buffer_cspc != None:
+                        if self.__nChannels > 1:
                             self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
                     #data_spc = numpy.sum(self.__buffer_spc,axis=0)
                     data_spc = numpy.nansum(self.__buffer_spc,axis=0)
-                    if self.__buffer_cspc != None:
+                    if self.__nChannels > 1:
                         #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
                         data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
                     else:

schainpy/model/proc/jroproc_voltage.py +403 -105

This diff has been collapsed as it changes many lines, (508 lines changed) Show them Hide them
	@@ -414,6 +414,14 class printAttribute(Operation):
	414	log.log(getattr(dataOut, attr), attr)	414	log.log(getattr(dataOut, attr), attr)
	415		415
	416	class cleanHeightsInterf(Operation):	416	class cleanHeightsInterf(Operation):
			417	__slots__ =('heights_indx', 'repeats', 'step', 'factor', 'idate', 'idxs','config','wMask')
			418	def __init__(self):
			419	self.repeats = 0
			420	self.factor=1
			421	self.wMask = None
			422	self.config = False
			423	self.idxs = None
			424	self.heights_indx = None
	417		425
	418	def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None):	426	def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None):
	419		427
	@@ -425,31 +433,53 class cleanHeightsInterf(Operation):
	425		433
	426	if currentTime < startTime or currentTime > endTime:	434	if currentTime < startTime or currentTime > endTime:
	427	return dataOut	435	return dataOut
			436	if not self.config:
			437
			438	#print(wMask)
			439	heights = [float(hei) for hei in heightsList]
			440	for r in range(repeats):
			441	heights += [ (h+(step*(r+1))) for h in heights]
			442	#print(heights)
			443	heiList = dataOut.heightList
			444	self.heights_indx = [getHei_index(h,h,heiList)[0] for h in heights]
			445
			446	self.wMask = numpy.asarray(factor)
			447	self.wMask = numpy.tile(self.wMask,(repeats+2))
			448	self.config = True
			449
			450	"""
			451	getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None)
			452	"""
			453	#print(self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref)))
	428		454
	429	wMask = numpy.asarray(factor)
	430	wMask = numpy.tile(wMask,(repeats+2))
	431	#print(wMask)
	432	heights = [float(hei) for hei in heightsList]
	433	for r in range(repeats):
	434	heights += [ (h+(step*(r+1))) for h in heights]
	435	#print(heights)
	436	heiList = dataOut.heightList
	437	heights_indx = [getHei_index(h,h,heiList)[0] for h in heights]
	438		455
	439	for ch in range(dataOut.data.shape[0]):	456	for ch in range(dataOut.data.shape[0]):
	440	i = 0	457	i = 0
	441	for hei in heights_indx:	458
	442	#print(dataOut.data[ch,hei])	459
			460	for hei in self.heights_indx:
			461	h = hei - 1
			462
			463
	443	if dataOut.data.ndim < 3:	464	if dataOut.data.ndim < 3:
	444	~~dataOut~~.~~data~~[ch,~~hei~~-1] = (dataOut.data[ch,hei-1])*~~wMask~~[i]	465	module = numpy.absolute(dataOut.data[ch,h])
			466	prev_h1 = numpy.absolute(dataOut.data[ch,h-1])
			467	dataOut.data[ch,h] = (dataOut.data[ch,h])/module * prev_h1
			468
			469	#dataOut.data[ch,hei-1] = (dataOut.data[ch,hei-1])*self.wMask[i]
	445	else:	470	else:
	446	~~dataOut~~.~~data~~[ch~~,:,hei~~-1] = (dataOut.data[ch,:,hei-1])*~~wMask~~[i]	471	module = numpy.absolute(dataOut.data[ch,:,h])
			472	prev_h1 = numpy.absolute(dataOut.data[ch,:,h-1])
			473	dataOut.data[ch,:,h] = (dataOut.data[ch,:,h])/module * prev_h1
			474	#dataOut.data[ch,:,hei-1] = (dataOut.data[ch,:,hei-1])*self.wMask[i]
	447	#print("done")	475	#print("done")
	448	i += 1	476	i += 1
	449		477
			478
	450	return dataOut	479	return dataOut
	451		480
	452		481
			482
	453	class interpolateHeights(Operation):	483	class interpolateHeights(Operation):
	454		484
	455	def run(self, dataOut, topLim, botLim):	485	def run(self, dataOut, topLim, botLim):
	@@ -1991,16 +2021,17 class removeProfileByFaradayHS(Operation):
	1991	outliers_IDs=[]	2021	outliers_IDs=[]
	1992	for c in range(nChannels):	2022	for c in range(nChannels):
	1993	for h in range(self.minHei_idx, self.maxHei_idx):	2023	for h in range(self.minHei_idx, self.maxHei_idx):
	1994	power = data[c,:,h] * numpy.conjugate(data[c,:,h])	2024	power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
	1995	power = power.real	2025	#power = power.real
	1996	#power = (numpy.abs(data[c,:,h].real))	2026	#power = (numpy.abs(data[c,:,h].real))
	1997	sortdata = numpy.sort(power, axis=None)	2027	sortdata = numpy.sort(power, axis=None)
	1998	sortID=power.argsort()	2028	sortID=power.argsort()
	1999	index = _noise.hildebrand_sekhon2(sortdata,self.navg) ~~#0.75-> buen valor~~	2029	index = _noise.hildebrand_sekhon2(sortdata,self.navg)
	2000		2030
	2001	indexes.append(index)	2031	indexes.append(index)
	2002	outliers_IDs=numpy.append(outliers_IDs,sortID[index:])	2032	outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
	2003	# print(outliers_IDs)	2033
			2034	# print(sortdata.min(), sortdata.max(), sortdata.mean())
	2004	# fig,ax = plt.subplots()	2035	# fig,ax = plt.subplots()
	2005	# #ax.set_title(str(k)+" "+str(j))	2036	# #ax.set_title(str(k)+" "+str(j))
	2006	# x=range(len(sortdata))	2037	# x=range(len(sortdata))
	@@ -2018,8 +2049,8 class removeProfileByFaradayHS(Operation):
	2018		2049
	2019		2050
	2020	#Agrupando el histograma de outliers,	2051	#Agrupando el histograma de outliers,
	2021	#my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)	2052	my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
	2022	my_bins = numpy.linspace(0,9600, 96, endpoint=False)	2053	#my_bins = numpy.linspace(0,1600, 96, endpoint=False)
	2023		2054
	2024	hist, bins = numpy.histogram(outs_lines,bins=my_bins)	2055	hist, bins = numpy.histogram(outs_lines,bins=my_bins)
	2025	hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier	2056	hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
	@@ -2040,22 +2071,104 class removeProfileByFaradayHS(Operation):
	2040		2071
	2041		2072
	2042		2073
	2043	# x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)	2074	x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
	2044	# fig, ax = plt.subplots(1,2,figsize=(8, 6))	2075	fig, ax = plt.subplots(1,2,figsize=(8, 6))
	2045	#	2076
	2046	# dat = data[0,:,:].real	2077	dat = data[0,:,:].real
	2047	# m = numpy.nanmean(dat)	2078	dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
	2048	~~# o~~ = numpy.nan~~std~~(dat)	2079	m = numpy.nanmean(dat)
	2049	# #print(m, o, x.shape, y.shape)	2080	o = numpy.nanstd(dat)
	2050	# c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2o), vmax = (m+2o))	2081	#print(m, o, x.shape, y.shape)
	2051	# ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')	2082	c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2o), vmax = (m+2o))
	2052	# fig.colorbar(c)	2083	ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
	2053	# ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')	2084	fig.colorbar(c)
	2054	# ax[1].hist(outs_lines,bins=my_bins)	2085	ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
	2055	# plt.show()	2086	ax[1].hist(outs_lines,bins=my_bins)
			2087	plt.show()
			2088
			2089
			2090	self.outliers_IDs_list = numpy.unique(outlier_loc_index)
			2091	print("outs list: ", self.outliers_IDs_list)
			2092	return data
			2093
			2094	def filterSatsProfiles2(self):
			2095	data = self.__buffer_data
			2096	#print(data.shape)
			2097	nChannels, profiles, heights = data.shape
			2098	indexes=numpy.zeros([], dtype=int)
			2099	outliers_IDs=[]
			2100	for c in range(nChannels):
			2101	noise_ref =10* numpy.log10((data[c,:,550:600] * numpy.conjugate(data[c,:,550:600])).real)
			2102	print("Noise ",noise_ref.mean())
			2103	for h in range(self.minHei_idx, self.maxHei_idx):
			2104	power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
			2105	#power = power.real
			2106	#power = (numpy.abs(data[c,:,h].real))
			2107	#sortdata = numpy.sort(power, axis=None)
			2108	#sortID=power.argsort()
			2109	#print(sortID)
			2110	th = 60 + 10
			2111	index = numpy.where(power > th )
			2112	if index[0].size > 10 and index[0].size < int(0.8*profiles):
			2113	indexes = numpy.append(indexes, index[0])
			2114	#print(index[0])
			2115	#print(index[0])
			2116
			2117	# fig,ax = plt.subplots()
			2118	# #ax.set_title(str(k)+" "+str(j))
			2119	# x=range(len(power))
			2120	# ax.scatter(x,power)
			2121	# #ax.axvline(index)
			2122	# plt.grid()
			2123	# plt.show()
			2124	#print(indexes)
			2125
			2126	#outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
			2127	#outliers_IDs = numpy.unique(outliers_IDs)
			2128
			2129	outs_lines = numpy.unique(indexes)
			2130	print("outliers Ids: ", outs_lines, outs_lines.shape)
			2131	#hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
			2132
			2133
			2134	#Agrupando el histograma de outliers,
			2135	my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
			2136	#my_bins = numpy.linspace(0,1600, 96, endpoint=False)
			2137
			2138	hist, bins = numpy.histogram(outs_lines,bins=my_bins)
			2139	hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
			2140	#print(hist_outliers_indexes[0])
			2141	bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
			2142	#print(bins_outliers_indexes)
			2143	outlier_loc_index = []
			2144
			2145
			2146
			2147	outlier_loc_index = [e for n in range(len(bins_outliers_indexes)-1) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin) ]
			2148
			2149	outlier_loc_index = numpy.asarray(outlier_loc_index)
			2150	outlier_loc_index = outlier_loc_index[~numpy.all(outlier_loc_index < 0)]
			2151
			2152	print("outliers final: ", outlier_loc_index)
			2153
			2154	x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
			2155	fig, ax = plt.subplots(1,2,figsize=(8, 6))
			2156
			2157	dat = data[0,:,:].real
			2158	dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
			2159	m = numpy.nanmean(dat)
			2160	o = numpy.nanstd(dat)
			2161	#print(m, o, x.shape, y.shape)
			2162	c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2o), vmax = (m+2o))
			2163	ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
			2164	fig.colorbar(c)
			2165	ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
			2166	ax[1].hist(outs_lines,bins=my_bins)
			2167	plt.show()
	2056		2168
	2057		2169
	2058	self.outliers_IDs_list = numpy.unique(outlier_loc_index)	2170	self.outliers_IDs_list = numpy.unique(outlier_loc_index)
			2171	print("outs list: ", self.outliers_IDs_list)
	2059	return data	2172	return data
	2060		2173
	2061	def cleanSpikesFFT2D(self):	2174	def cleanSpikesFFT2D(self):
	@@ -2588,7 +2701,7 class removeProfileByFaradayHS(Operation):
	2588	#print("apnd : ",data)	2701	#print("apnd : ",data)
	2589	#dataBlock = self.cleanOutliersByBlock()	2702	#dataBlock = self.cleanOutliersByBlock()
	2590	#dataBlock = self.cleanSpikesFFT2D()	2703	#dataBlock = self.cleanSpikesFFT2D()
	2591	dataBlock = self.filterSatsProfiles()	2704	dataBlock = self.filterSatsProfiles2()
	2592	self.__dataReady = True	2705	self.__dataReady = True
	2593		2706
	2594	return dataBlock	2707	return dataBlock
	@@ -2649,7 +2762,7 class removeProfileByFaradayHS(Operation):
	2649	#print("prof: ",self.init_prof)	2762	#print("prof: ",self.init_prof)
	2650	dataOut.flagNoData = False	2763	dataOut.flagNoData = False
	2651	if numpy.isin(self.n_prof_released, self.outliers_IDs_list):	2764	if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
	2652	#print("omitting: ", self.n_prof_released)	2765	print("omitting: ", self.n_prof_released)
	2653	dataOut.flagNoData = True	2766	dataOut.flagNoData = True
	2654	dataOut.ippSeconds = self._ipp	2767	dataOut.ippSeconds = self._ipp
	2655	dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp	2768	dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
	@@ -2664,7 +2777,7 class removeProfileByFaradayHS(Operation):
	2664		2777
	2665	self.init_prof = self.end_prof	2778	self.init_prof = self.end_prof
	2666	self.end_prof += self.lenProfileOut	2779	self.end_prof += self.lenProfileOut
	2667	# #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)	2780	#print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
	2668	self.n_prof_released += 1	2781	self.n_prof_released += 1
	2669		2782
	2670	if self.end_prof >= (self.n +self.lenProfileOut):	2783	if self.end_prof >= (self.n +self.lenProfileOut):
	@@ -2722,9 +2835,10 class removeProfileByFaradayHS(Operation):
	2722		2835
	2723	return dataOut	2836	return dataOut
	2724		2837
			2838
	2725	class RemoveProfileSats(Operation):	2839	class RemoveProfileSats(Operation):
	2726	'''	2840	'''
	2727	Omite los perfiles contaminados con señal de satelites,	2841	Omite los perfiles contaminados con señal de satélites, usando una altura de referencia
	2728	In: minHei = min_sat_range	2842	In: minHei = min_sat_range
	2729	max_sat_range	2843	max_sat_range
	2730	min_hei_ref	2844	min_hei_ref
	@@ -2734,106 +2848,290 class RemoveProfileSats(Operation):
	2734	profile clean	2848	profile clean
	2735	'''	2849	'''
	2736		2850
	2737	isConfig = False
	2738	min_sats = 0
	2739	max_sats = 999999999
	2740	min_ref= 0
	2741	max_ref= 9999999999
	2742	needReshape = False
	2743	count = 0
	2744	thdB = 0
	2745	byRanges = False
	2746	min_sats = None
	2747	max_sats = None
	2748	noise = 0
	2749		2851
			2852	__buffer_data = []
			2853	__buffer_times = []
			2854
			2855	buffer = None
			2856
			2857	outliers_IDs_list = []
			2858
			2859
			2860	__slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
			2861	'first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
			2862	'__count_exec','__initime','__dataReady','__ipp', 'minRef', 'maxRef', 'thdB')
	2750	def __init__(self, **kwargs):	2863	def __init__(self, **kwargs):
	2751		2864
	2752	Operation.__init__(self, **kwargs)	2865	Operation.__init__(self, **kwargs)
	2753	self.isConfig = False	2866	self.isConfig = False
	2754		2867
			2868	def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3,
			2869	minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
			2870
			2871	if n == None and timeInterval == None:
			2872	raise ValueError("nprofiles or timeInterval should be specified ...")
	2755		2873
	2756	def setup(self, dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList):	2874	if n != None:
			2875	self.n = n
	2757		2876
	2758	if rangeHeiList!=None:	2877	self.navg = navg
	2759	self.byRanges = True	2878	self.profileMargin = profileMargin
	2760	else:	2879	self.thHistOutlier = thHistOutlier
	2761	if minHei==None or maxHei==None :	2880	self.__profIndex = 0
	2762	raise ValueError("Parameters heights are required")	2881	self.buffer = None
	2763	if minRef==None or maxRef==None:	2882	self._ipp = dataOut.ippSeconds
	2764	raise ValueError("Parameters heights are required")	2883	self.n_prof_released = 0
	2765		2884	self.heightList = dataOut.heightList
	2766	if self.byRanges:	2885	self.init_prof = 0
	2767	self.~~min_sats~~ = []	2886	self.end_prof = 0
	2768	self.max_sats = []	2887	self.__count_exec = 0
	2769	for min,max in rangeHeiList:	2888	self.__profIndex = 0
	2770	a,b = getHei_index(min, max, dataOut.heightList)	2889	self.first_utcBlock = None
	2771	self.min_sats.append(a)	2890	#self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
	2772	self.max_sats.append(b)	2891	minHei = minHei
	2773	else:	2892	maxHei = maxHei
	2774	self.min_sats, self.max_sats = getHei_index(minHei, maxHei, dataOut.heightList)	2893	if minHei==None :
			2894	minHei = dataOut.heightList[0]
			2895	if maxHei==None :
			2896	maxHei = dataOut.heightList[-1]
			2897	self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
	2775	self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)	2898	self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
	2776	self.th = th	2899	self.nChannels = dataOut.nChannels
			2900	self.nHeights = dataOut.nHeights
			2901	self.test_counter = 0
	2777	self.thdB = thdB	2902	self.thdB = thdB
	2778	self.isConfig = True
	2779		2903
			2904	def filterSatsProfiles(self):
			2905	data = self.__buffer_data
			2906	#print(data.shape)
			2907	nChannels, profiles, heights = data.shape
			2908	indexes=numpy.zeros([], dtype=int)
			2909	outliers_IDs=[]
			2910	for c in range(nChannels):
			2911	#print(self.min_ref,self.max_ref)
			2912	noise_ref = 10* numpy.log10((data[c,:,self.min_ref:self.max_ref] * numpy.conjugate(data[c,:,self.min_ref:self.max_ref])).real)
			2913	#print("Noise ",numpy.percentile(noise_ref,95))
			2914	p95 = numpy.percentile(noise_ref,95)
			2915	noise_ref = noise_ref.mean()
			2916	#print("Noise ",noise_ref
	2780		2917
	2781	def compareRanges(self,data, minHei,maxHei):
	2782		2918
	2783	# ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref])	2919	for h in range(self.minHei_idx, self.maxHei_idx):
	2784	# p_ref = 10*numpy.log10(ref.real)	2920	power = 10* numpy.log10((data[c,:,h] * numpy.conjugate(data[c,:,h])).real)
	2785	# m_ref = numpy.mean(p_ref)	2921	#th = noise_ref + self.thdB
			2922	th = noise_ref + 1.5*(p95-noise_ref)
			2923	index = numpy.where(power > th )
			2924	if index[0].size > 10 and index[0].size < int(self.navg*profiles):
			2925	indexes = numpy.append(indexes, index[0])
			2926	#print(index[0])
			2927	#print(index[0])
			2928
			2929	# fig,ax = plt.subplots()
			2930	# #ax.set_title(str(k)+" "+str(j))
			2931	# x=range(len(power))
			2932	# ax.scatter(x,power)
			2933	# #ax.axvline(index)
			2934	# plt.grid()
			2935	# plt.show()
			2936	#print(indexes)
	2786		2937
	2787	m_ref = self.noise	2938	#outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
			2939	#outliers_IDs = numpy.unique(outliers_IDs)
	2788		2940
	2789	sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])	2941	outs_lines = numpy.unique(indexes)
	2790	p_sats = 10*numpy.log10(sats.real)
	2791	m_sats = numpy.mean(p_sats)
	2792		2942
	2793	if m_sats > (m_ref + self.th): #and (m_sats > self.thdB):
	2794	#print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref))
	2795	#print("Removing profiles...")
	2796	return False
	2797		2943
	2798	return True	2944	#Agrupando el histograma de outliers,
			2945	my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=True)
	2799		2946
	2800	def isProfileClean(self, data):
	2801	'''
	2802	Analiza solo 1 canal, y descarta todos...
	2803	'''
	2804		2947
	2805	clean = True	2948	hist, bins = numpy.histogram(outs_lines,bins=my_bins)
			2949	hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
			2950	hist_outliers_indexes = hist_outliers_indexes[0]
			2951	if len(hist_outliers_indexes>0):
			2952	hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1)
			2953	#print(hist_outliers_indexes)
			2954	#print(bins, hist_outliers_indexes)
			2955	bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] #
			2956
			2957	outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n]+ profiles//100 + self.profileMargin) ]
			2958	outlier_loc_index = numpy.asarray(outlier_loc_index)
			2959
			2960	#print("outliers Ids: ", outlier_loc_index, outlier_loc_index.shape)
			2961	outlier_loc_index = outlier_loc_index[ (outlier_loc_index >= 0) & (outlier_loc_index<profiles)]
			2962	#print("outliers final: ", outlier_loc_index)
			2963
			2964	# from matplotlib import pyplot as plt
			2965	# x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
			2966	# fig, ax = plt.subplots(1,2,figsize=(8, 6))
			2967	# dat = data[0,:,:].real
			2968	# dat = 10* numpy.log10((data[0,:,:] * numpy.conjugate(data[0,:,:])).real)
			2969	# m = numpy.nanmean(dat)
			2970	# o = numpy.nanstd(dat)
			2971	# #print(m, o, x.shape, y.shape)
			2972	# c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2o), vmax = (m+2o))
			2973	# ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
			2974	# fig.colorbar(c)
			2975	# ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
			2976	# ax[1].hist(outs_lines,bins=my_bins)
			2977	# plt.show()
	2806		2978
	2807	if self.byRanges:
	2808		2979
	2809	for n in range(len(self.min_sats)):	2980	self.outliers_IDs_list = outlier_loc_index
	2810	c = self.compareRanges(data,self.min_sats[n],self.max_sats[n])	2981	#print("outs list: ", self.outliers_IDs_list)
	2811	clean = clean and c	2982	return data
			2983
			2984
			2985
			2986	def fillBuffer(self, data, datatime):
			2987
			2988	if self.__profIndex == 0:
			2989	self.__buffer_data = data.copy()
			2990
	2812	else:	2991	else:
			2992	self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
			2993	self.__profIndex += 1
			2994	self.__buffer_times.append(datatime)
	2813		2995
	2814	clean = (self.compareRanges(data, self.min_sats,self.max_sats))	2996	def getData(self, data, datatime=None):
	2815		2997
	2816	return clean	2998	if self.__profIndex == 0:
			2999	self.__initime = datatime
	2817		3000
	2818		3001
			3002	self.__dataReady = False
			3003
			3004	self.fillBuffer(data, datatime)
			3005	dataBlock = None
			3006
			3007	if self.__profIndex == self.n:
			3008	#print("apnd : ",data)
			3009	dataBlock = self.filterSatsProfiles()
			3010	self.__dataReady = True
			3011
			3012	return dataBlock
			3013
			3014	if dataBlock is None:
			3015	return None, None
	2819		3016
	2820	def run(self, dataOut, minHei=None, maxHei=None, minRef=None, maxRef=None, th=5, thdB=65, rangeHeiList=None):	3017
	2821	dataOut.flagNoData = True	3018
			3019	return dataBlock
			3020
			3021	def releaseBlock(self):
			3022
			3023	if self.n % self.lenProfileOut != 0:
			3024	raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
			3025	return None
			3026
			3027	data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
			3028
			3029	self.init_prof = self.end_prof
			3030	self.end_prof += self.lenProfileOut
			3031	#print("data release shape: ",dataOut.data.shape, self.end_prof)
			3032	self.n_prof_released += 1
			3033
			3034	return data
			3035
			3036	def run(self, dataOut, n=None, navg=0.75, nProfilesOut=1, profile_margin=50,
			3037	th_hist_outlier=3,minHei=None, maxHei=None, minRef=None, maxRef=None, thdB=10):
	2822		3038
	2823	if not self.isConfig:	3039	if not self.isConfig:
	2824	self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)	3040	#print("init p idx: ", dataOut.profileIndex )
			3041	self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,thHistOutlier=th_hist_outlier,
			3042	minHei=minHei, maxHei=maxHei, minRef=minRef, maxRef=maxRef, thdB=thdB)
	2825	self.isConfig = True	3043	self.isConfig = True
	2826	#print(self.min_sats,self.max_sats)
	2827	if dataOut.flagDataAsBlock:
	2828	raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
	2829		3044
	2830	else:	3045	dataBlock = None
	2831	self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref))	3046
	2832	if not self.isProfileClean(dataOut.data):	3047	if not dataOut.buffer_empty: #hay datos acumulados
	2833	return dataOut	3048
	2834	#dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)	3049	if self.init_prof == 0:
	2835	#self.count += 1	3050	self.n_prof_released = 0
			3051	self.lenProfileOut = nProfilesOut
			3052	dataOut.flagNoData = False
			3053	#print("tp 2 ",dataOut.data.shape)
			3054
			3055	self.init_prof = 0
			3056	self.end_prof = self.lenProfileOut
			3057
			3058	dataOut.nProfiles = self.lenProfileOut
			3059	if nProfilesOut == 1:
			3060	dataOut.flagDataAsBlock = False
			3061	else:
			3062	dataOut.flagDataAsBlock = True
			3063	#print("prof: ",self.init_prof)
			3064	dataOut.flagNoData = False
			3065	if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
			3066	#print("omitting: ", self.n_prof_released)
			3067	dataOut.flagNoData = True
			3068	dataOut.ippSeconds = self._ipp
			3069	dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
			3070	# print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
			3071	#dataOut.data = self.releaseBlock()
			3072	#########################################################3
			3073	if self.n % self.lenProfileOut != 0:
			3074	raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
			3075	return None
			3076
			3077	dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
			3078
			3079	self.init_prof = self.end_prof
			3080	self.end_prof += self.lenProfileOut
			3081	#print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
			3082	self.n_prof_released += 1
			3083
			3084	if self.end_prof >= (self.n +self.lenProfileOut):
			3085
			3086	self.init_prof = 0
			3087	self.__profIndex = 0
			3088	self.buffer = None
			3089	dataOut.buffer_empty = True
			3090	self.outliers_IDs_list = []
			3091	self.n_prof_released = 0
			3092	dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
			3093	#print("cleaning...", dataOut.buffer_empty)
			3094	dataOut.profileIndex = 0 #self.lenProfileOut
			3095	####################################################################
			3096	return dataOut
			3097
	2836		3098
	2837	dataOut.flagNoData = False	3099	#print("tp 223 ",dataOut.data.shape)
			3100	dataOut.flagNoData = True
			3101
			3102
			3103
			3104	try:
			3105	#dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
			3106	dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
			3107	self.__count_exec +=1
			3108	except Exception as e:
			3109	print("Error getting profiles data",self.__count_exec )
			3110	print(e)
			3111	sys.exit()
			3112
			3113	if self.__dataReady:
			3114	#print("omitting: ", len(self.outliers_IDs_list))
			3115	self.__count_exec = 0
			3116	#dataOut.data =
			3117	#self.buffer = numpy.flip(dataBlock, axis=1)
			3118	self.buffer = dataBlock
			3119	self.first_utcBlock = self.__initime
			3120	dataOut.utctime = self.__initime
			3121	dataOut.nProfiles = self.__profIndex
			3122	#dataOut.flagNoData = False
			3123	self.init_prof = 0
			3124	self.__profIndex = 0
			3125	self.__initime = None
			3126	dataBlock = None
			3127	self.__buffer_times = []
			3128	dataOut.error = False
			3129	dataOut.useInputBuffer = True
			3130	dataOut.buffer_empty = False
			3131	#print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
			3132
			3133
			3134
			3135	#print(self.__count_exec)
	2838		3136
	2839	return dataOut	3137	return dataOut

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