schain Commit - r1546:d5ba7b76d333 · Jicamarca Repository

isr procs ready

joabAM -

r1546:d5ba7b76d333

parent child

Context file:

r1546:d5ba7b76d333

Collapse all files

schainpy/model/proc/jroproc_spectra.py +30 -35

                     self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
                     self.isConfig = True
+                    self.offset = 1
                     if offset!=None:
                         self.offset = 10**(offset/10)
                     #print("config getNoiseB Done")
                     #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
                     #print(self.dataOut.flagNoData)
-                    print("getNoise Done", noise)
+                    #print("getNoise Done", noise)
                     return self.dataOut
                 def getNoiseByMean(self,data):
                     buffer=None
                     buffer1=None
                     buffer_cspc=None
+                    #print("aes: ", self.__buffer_cspc)
                     self.__buffer_spc=numpy.array(self.__buffer_spc)
-                    try:
+                    if self.__buffer_cspc != None:
                         self.__buffer_cspc=numpy.array(self.__buffer_cspc)
-                    except :
-                        #print("No cpsc",e)
-                        pass
                     #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
                     freq_dc = int(self.__buffer_spc.shape[2] / 2)
                     self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
                     for k in range(self.minHei_ind,self.maxHei_ind):
-                        try:
+                        if self.__buffer_cspc != None:
                             buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
-                        except:
-                            #print("No cpsc",e)
-                            pass
                         outliers_IDs_cspc=[]
                         cspc_outliers_exist=False
                         for i in range(self.nChannels):#dataOut.nChannels):
                                 self.dataOutliers[i,k] = len(outliers_IDs)
                             self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
-                            outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
+                            if self.__buffer_cspc != None:
+                                outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
-                        outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
+                        if self.__buffer_cspc != None:
-                        if cspc_outliers_exist :
+                            outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
+                        if cspc_outliers_exist:
                             lt=outliers_IDs_cspc
                                 #buffer_cspc[p,:]=avg
                                 buffer_cspc[p,:] = numpy.NaN
-                        try:
+                        if self.__buffer_cspc != None:
                             self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
-                        except:
-                            #print("No cpsc",e)
-                            pass
-                        #else:
-                            #break
                     #data_spc = numpy.sum(self.__buffer_spc,axis=0)
                     data_spc = numpy.nansum(self.__buffer_spc,axis=0)
-                    try:
+                    if self.__buffer_cspc != None:
                         #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
                         data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
-                    except:
+                    else:
-                        #print("No cpsc",e)
+                        data_cspc = None
-                        pass
                     data_dc = self.__buffer_dc
                     #(CH, HEIGH)
                     self.maxProfilesInt = self.__profIndex
                     self.__buffer_cspc = []
                     self.__buffer_dc = 0
                     self.__profIndex = 0
+                    #print("cleaned ",data_cspc)
                     return data_spc, data_cspc, data_dc, n
                 def byProfiles(self, *args):
                     if not self.__dataReady:
                         return None, None, None, None
+                    #print("integrate", avgdata_cspc)
                     return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
                 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
                     if n == 1:
                         return self.dataOut
+                    #print("nchannels", self.dataOut.nChannels)
                     if self.dataOut.nChannels == 1:
                         self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
-                    #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc)
+                    #print("IN spc:", self.dataOut.data_spc.shape, self.dataOut.data_cspc)
                     if not self.isConfig:
                         self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
                         self.isConfig = True
                             if self.nChannels == 1:
                                 #print("f int", avgdata_spc.shape)
                                 self.dataOut.data_spc = avgdata_spc
-                                self.dataOut.data_cspc = avgdata_spc
+                                self.dataOut.data_cspc = None
                             else:
                                 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
                                 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
                         #print(self.dataOut.max_nIncohInt)
                         self.dataOut.utctime = avgdatatime
                         self.dataOut.flagNoData = False
-                        #print("Faraday Integration DONE...")
+                        #print("Faraday Integration DONE...", self.dataOut.data_cspc)
                     #print(self.dataOut.flagNoData)
                     return self.dataOut
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     jnoise = self.dataOut.getNoise()
-                    num_incoh = self.dataOut.nIncohInt
+                    #num_incoh = self.dataOut.nIncohInt
+                    num_incoh = self.dataOut.max_nIncohInt
+                    #print("spc: ", jspectra.shape, jcspectra)
                     num_channel = jspectra.shape[0]
                     num_prof = jspectra.shape[1]
                     num_hei = jspectra.shape[2]
                         count_hei = int(num_hei / 2)
                         hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
                         hei_interf = numpy.asarray(hei_interf)[0]
+                        #print(hei_interf)
                     # nhei_interf
                     if (nhei_interf == None):
                         nhei_interf = 5
                         jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
                         jspc_interf = jspc_interf.transpose()
                         # Calculando el espectro de interferencia promedio
-                        noiseid = numpy.where(
+                        noiseid = numpy.where(jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
-                            jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
-                        interfid = numpy.where(
+                        interfid = numpy.where(jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
-                            jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
                         interfid = interfid[0]
                         cinterfid = interfid.size

schainpy/model/proc/jroproc_voltage.py +109 -13

@@ -2023,6 +2023,7 class removeProfileByFaradayHS(Operation):
2023	def cleanOutliersByBlock(self):	2023	def cleanOutliersByBlock(self):
2024	import matplotlib.pyplot as plt	2024	import matplotlib.pyplot as plt
2025	import datetime	2025	import datetime
		2026	import cv2
2026	#print(self.__buffer_data[0].shape)	2027	#print(self.__buffer_data[0].shape)
2027	data = self.__buffer_data#.copy()	2028	data = self.__buffer_data#.copy()
2028	print("cleaning shape inpt: ",data.shape)	2029	print("cleaning shape inpt: ",data.shape)
@@ -2034,59 +2035,154 class removeProfileByFaradayHS(Operation):
2034	(nch,nsamples, nh) = spectrum.shape	2035	(nch,nsamples, nh) = spectrum.shape
2035	data2 = None	2036	data2 = None
2036	#print(data.shape)	2037	#print(data.shape)
		2038	cleanedBlock = None
2037	spectrum2 = spectrum.copy()	2039	spectrum2 = spectrum.copy()
2038	for ch in range(nch):	2040	for ch in range(nch):
2039	dh = self.__dh	2041	dh = self.__dh
2040	dt1 = (dh*1.0e-6)/(0.15)	2042	dt1 = (dh*1.0e-6)/(0.15)
2041	dt2 = self.__buffer_times[1]-self.__buffer_times[0]	2043	dt2 = self.__buffer_times[1]-self.__buffer_times[0]
		2044
2042	freqv = numpy.fft.fftfreq(nh, d=dt1)	2045	freqv = numpy.fft.fftfreq(nh, d=dt1)
2043	freqh = numpy.fft.fftfreq(self.n, d=dt2)	2046	freqh = numpy.fft.fftfreq(self.n, d=dt2)
2044	#print("spc loop: ")	2047
		2048	# freqv = numpy.fft.fftshift(numpy.fft.fftfreq(nh, d=dt1))
		2049	# freqh = numpy.fft.fftshift(numpy.fft.fftfreq(self.n, d=dt2))
		2050	#
		2051
		2052
		2053
2045		2054
2046		2055
2047		2056
2048	x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
2049	z = numpy.abs(spectrum[ch,:,:])	2057	z = numpy.abs(spectrum[ch,:,:])
2050	phase = numpy.angle(spectrum[ch,:,:])	2058	phase = numpy.angle(spectrum[ch,:,:])
		2059
		2060	print("shape z: ", z.shape)
		2061
		2062
2051	# Find all peaks higher than the 98th percentile	2063	# Find all peaks higher than the 98th percentile
2052	peaks = z < numpy.percentile(z, 99)	2064	peaks = z < numpy.percentile(z, 99)
2053	#print(peaks)	2065	#print(peaks)
2054	# Set those peak coefficients to zero	2066	# Set those peak coefficients to zero
2055	spectrum2 = spectrum2 * peaks.astype(int)	2067	#spectrum2 = spectrum * peaks.astype(int)
2056	data2 = numpy.fft.ifft2(spectrum2,axes=(0,2))	2068	#data2 = numpy.fft.ifft2(spectrum2,axes=(0,2))
2057		2069
2058	dat = numpy.log10(z.T)	2070	dat = numpy.log10(z.T)
2059	pdat = numpy.log10(phase.T)	2071	pdat = numpy.log10(phase.T)
2060	dat2 = numpy.log10(numpy.abs(spectrum2.T))	2072
		2073
		2074	min, max = dat.min(), dat.max()
		2075	#img = ((dat-min)*255/(max-min)).astype(numpy.uint8)
		2076	img = ((dat-min)*200/(max-min)).astype(numpy.uint8)
		2077	#img = 255 - ((dat-min)*255/max).astype(numpy.uint8)
		2078
		2079	#print(img.shape)
		2080	# cv2.imshow('image', numpy.fft.fftshift(img))
		2081	# cv2.waitKey(0)
		2082
		2083	spikes = numpy.where(img > 230, True, False)
		2084	#print("spikes: ", spikes.sum())
		2085
		2086	#img = cv2.medianBlur(img,3)
		2087
		2088	#cv2.imshow('image', cv2.resize(img.astype('float32'), (600, 600)))
		2089	kernel3 = numpy.zeros((3,15))
		2090	kernel3[0, :] = -1
		2091	kernel3[1, :] = 2
		2092	kernel3[2, :] = -1
		2093	sharp_img = cv2.filter2D(src=img, ddepth=-1, kernel=kernel3)
		2094
		2095
		2096	# cv2.imshow('sharp image', numpy.fft.fftshift(sharp_img))
		2097	# cv2.waitKey(0)
		2098
		2099	ret, thresh = cv2.threshold(sharp_img, 127, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
		2100	# #thresh = cv2.adaptiveThreshold(sharp_img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY_INV,21,0)
		2101	# cv2.imshow('binary img', numpy.fft.fftshift(thresh))
		2102	# cv2.waitKey(0)
		2103	# #Remove horizontal
		2104	horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (10,1))
		2105	detected_lines = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=8)
		2106	cnts, h = cv2.findContours(detected_lines, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
		2107	# print(cnts)
		2108	# cv2.imshow('lines ', numpy.fft.fftshift(detected_lines))
		2109	# cv2.waitKey(0)
		2110	#cnts = cnts[0] if len(cnts) == 2 else cnts[1]
		2111	y_line_index = numpy.asarray([ [c[0][0][0],c[1][0][0], c[0][0][1]] for c in cnts ])
		2112
		2113	# print(y_line_index)
		2114	u, index, counts = numpy.unique(y_line_index[:,2], return_index=True, return_counts=True)
		2115	print(u, index, counts)
		2116	readed_lines = []
		2117	for x1, x2, y, YY, n in zip(y_line_index, u, counts):
		2118	if y in readed_lines:
		2119	pass
		2120	if n > 1:
		2121	len_line = self.n
		2122	else:
		2123	len_line = len(z[x1:x2,y])
		2124	#print(y, nh, self.n)
		2125	if y != (nh-1):
		2126	list = [ ((z[n, y-1] + z[n,y+1])/2) for n in range(len_line)]
		2127	else:
		2128	list = [ ((z[n, y-1] + z[n,0])/2) for n in range(len_line)]
		2129	if n > 1:
		2130	z[:,y] = numpy.asarray(list)
		2131	else:
		2132	z[x1:x2,y] = numpy.asarray(list)
		2133
		2134	readed_lines.append(y)
		2135
		2136	dat2 = numpy.log10(z.T)
		2137	img = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
		2138
		2139	for c in cnts:
		2140	#print(c)
		2141	cv2.drawContours(img, [c], -1, (255,255,255), 2)
		2142	#print("C: ",c[0][1])
		2143
		2144	spcCleaned = z * numpy.exp(1j*phase)
		2145	#print(spcCleaned)
		2146
		2147
		2148	cv2.imshow('image2', numpy.fft.fftshift(img))
		2149	cv2.waitKey(0)
		2150	cv2.destroyAllWindows()
2061		2151
2062	m = numpy.mean(dat)	2152	m = numpy.mean(dat)
2063	o = numpy.std(dat)	2153	o = numpy.std(dat)
		2154
2064	fig, ax = plt.subplots(1,2,figsize=(12, 6))	2155	fig, ax = plt.subplots(1,2,figsize=(12, 6))
		2156	#X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
		2157	X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv))
2065		2158
2066	colormap = 'jet'	2159	colormap = 'jet'
2067	#c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2o)/2, vmax = (m+2o))	2160	#c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2o)/2, vmax = (m+2o))
2068	c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = 4.2, vmax = 5.0)	2161	c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 4.2, vmax = 5.0)
2069	fig.colorbar(c, ax=ax[0])	2162	fig.colorbar(c, ax=ax[0])
2070	#print("aqui estoy", dat.shape)	2163
		2164
2071	#c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2o), vmax = (m+2o))	2165	#c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2o), vmax = (m+2o))
2072	date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')	2166	#date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
2073	#strftime('%Y-%m-%d %H:%M:%S')	2167	#strftime('%Y-%m-%d %H:%M:%S')
2074	#ax[0].set_title('Spectrum magnitude '+date_time)	2168	#ax[0].set_title('Spectrum magnitude '+date_time)
2075	fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)	2169	#fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
2076	#print("aqui estoy2",dat2[:,:,0].shape)	2170	#print("aqui estoy2",dat2[:,:,0].shape)
2077	c = ax[1].pcolormesh(x, y, pdat, cmap =colormap, vmin = -~~0.0~~, vmax = 0.5)	2171	c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = 4.2, vmax = 5.0)
2078	#c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2o)/2, vmax = (m+2o)-1)	2172	#c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2o)/2, vmax = (m+2o)-1)
2079	#print("aqui estoy3")	2173	#print("aqui estoy3")
2080	fig.colorbar(c, ax=ax[1])	2174	fig.colorbar(c, ax=ax[1])
2081	plt.show()	2175	plt.show()
2082		2176
		2177	spectrum[ch,:,:] = spcCleaned
		2178
2083	#print(data2.shape)	2179	#print(data2.shape)
2084		2180
2085	#data = data2
2086		2181
2087	#cleanBlock = numpy.fft.ifft2(data, axes=(0,2)).reshape()
2088		2182
2089	#print("cleanOutliersByBlock Done")	2183	data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(0,2))
		2184
		2185	#print("cleanOutliersByBlock Done", data.shape)
2090		2186
2091	return data	2187	return data
2092		2188

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