##// END OF EJS Templates
schain
RSS
proc schain amisr isr1
joabAM -
r1547:f536cbdd0d0e
parent child
Show More
Side by Side Unified
Context file:
r1547:f536cbdd0d0e Loading diff...:
Show file before | Show file after | Hide comments
@@ -60,7 +60,7 class SpectraPlot(Plot):
60 data['rti'] = dataOut.getPower()
60 data['rti'] = dataOut.getPower()
61 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
61 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
62 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
62 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
63 data['noise'] = noise[0]
63 data['noise'] = noise
64
64
65 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
65 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
66 if self.CODE == 'spc_moments':
66 if self.CODE == 'spc_moments':
@@ -91,7 +91,7 class SpectraPlot(Plot):
91 z = data['spc']
91 z = data['spc']
92 #print(z.shape, x.shape, y.shape)
92 #print(z.shape, x.shape, y.shape)
93 for n, ax in enumerate(self.axes):
93 for n, ax in enumerate(self.axes):
94 noise = self.data['noise'][n]
94 noise = self.data['noise'][n][0]
95 #print(noise)
95 #print(noise)
96 if self.CODE == 'spc_moments':
96 if self.CODE == 'spc_moments':
97 mean = data['moments'][n, 1]
97 mean = data['moments'][n, 1]
@@ -275,7 +275,7 class RTIPlot(Plot):
275 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
275 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
276 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
276 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
277 data['noise'] = noise
277 data['noise'] = noise
278
278
279 return data, meta
279 return data, meta
280
280
281 def plot(self):
281 def plot(self):
@@ -533,9 +533,10 class SpectraCutPlot(Plot):
533 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
533 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
534 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
534 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
535
535
536 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - n0
536 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
537 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
537
538
538 data['spc'] = spc
539 data['spc'] = spc - noise
539 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
540 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
540
541
541 return data, meta
542 return data, meta
@@ -827,7 +828,7 class NoiselessSpectraPlot(Plot):
827 '''
828 '''
828
829
829 CODE = 'noiseless_spc'
830 CODE = 'noiseless_spc'
830 colormap = 'nipy_spectral'
831 colormap = 'jet'
831 plot_type = 'pcolor'
832 plot_type = 'pcolor'
832 buffering = False
833 buffering = False
833 channelList = []
834 channelList = []
@@ -866,8 +867,12 class NoiselessSpectraPlot(Plot):
866
867
867 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
868 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
868
869
869 data['spc'] = spc - n0
870 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
870 data['rti'] = dataOut.getPower() - n0
871 data['rti'] = dataOut.getPower() - noise
872
873 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
874 data['spc'] = spc - noise
875
871
876
872 # data['noise'] = noise
877 # data['noise'] = noise
873 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
878 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
@@ -963,11 +968,13 class NoiselessRTIPlot(Plot):
963
968
964
969
965 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
970 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
971 #print("Norm: ", norm)
972 #print(dataOut.nProfiles , dataOut.max_nIncohInt ,dataOut.nCohInt, dataOut.windowOfFilter)
966 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
973 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
967
974
968 data['noise'] = n0[0]
975 data['noise'] = n0
969
976 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
970 data['noiseless_rti'] = dataOut.getPower() - n0
977 data['noiseless_rti'] = dataOut.getPower() - noise
971
978
972 return data, meta
979 return data, meta
973
980
Show file before | Show file after | Hide comments
@@ -128,7 +128,7 class ParametersProc(ProcessingUnit):
128 #---------------------- Spectra Data ---------------------------
128 #---------------------- Spectra Data ---------------------------
129
129
130 if self.dataIn.type == "Spectra":
130 if self.dataIn.type == "Spectra":
131
131
132 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
132 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
133 self.dataOut.data_spc = self.dataIn.data_spc
133 self.dataOut.data_spc = self.dataIn.data_spc
134 self.dataOut.data_cspc = self.dataIn.data_cspc
134 self.dataOut.data_cspc = self.dataIn.data_cspc
Show file before | Show file after | Hide comments
@@ -22,6 +22,7 from schainpy.model.data import _noise
22 from schainpy.utils import log
22 from schainpy.utils import log
23 import matplotlib.pyplot as plt
23 import matplotlib.pyplot as plt
24 #from scipy.optimize import curve_fit
24 #from scipy.optimize import curve_fit
25 from schainpy.model.io.utils import getHei_index
25
26
26 class SpectraProc(ProcessingUnit):
27 class SpectraProc(ProcessingUnit):
27
28
@@ -73,6 +74,7 class SpectraProc(ProcessingUnit):
73
74
74
75
75 def __getFft(self):
76 def __getFft(self):
77 # print("fft donw")
76 """
78 """
77 Convierte valores de Voltaje a Spectra
79 Convierte valores de Voltaje a Spectra
78
80
@@ -124,7 +126,7 class SpectraProc(ProcessingUnit):
124 self.dataOut.blockSize = blocksize
126 self.dataOut.blockSize = blocksize
125 self.dataOut.flagShiftFFT = False
127 self.dataOut.flagShiftFFT = False
126
128
127 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
129 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False, zeroPad=False):
128 #print("run spc proc")
130 #print("run spc proc")
129 try:
131 try:
130 type = self.dataIn.type.decode("utf-8")
132 type = self.dataIn.type.decode("utf-8")
@@ -167,15 +169,21 class SpectraProc(ProcessingUnit):
167 self.dataOut.nFFTPoints = nFFTPoints
169 self.dataOut.nFFTPoints = nFFTPoints
168 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
170 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
169 if self.buffer is None:
171 if self.buffer is None:
170 self.buffer = numpy.zeros((self.dataIn.nChannels,
172 if not zeroPad:
173 self.buffer = numpy.zeros((self.dataIn.nChannels,
171 nProfiles,
174 nProfiles,
172 self.dataIn.nHeights),
175 self.dataIn.nHeights),
173 dtype='complex')
176 dtype='complex')
177 else:
178 self.buffer = numpy.zeros((self.dataIn.nChannels,
179 nFFTPoints,
180 self.dataIn.nHeights),
181 dtype='complex')
174
182
175 if self.dataIn.flagDataAsBlock:
183 if self.dataIn.flagDataAsBlock:
176 nVoltProfiles = self.dataIn.data.shape[1]
184 nVoltProfiles = self.dataIn.data.shape[1]
177
185
178 if nVoltProfiles == nProfiles:
186 if nVoltProfiles == nProfiles or zeroPad:
179 self.buffer = self.dataIn.data.copy()
187 self.buffer = self.dataIn.data.copy()
180 self.profIndex = nVoltProfiles
188 self.profIndex = nVoltProfiles
181
189
@@ -201,7 +209,7 class SpectraProc(ProcessingUnit):
201 if self.firstdatatime == None:
209 if self.firstdatatime == None:
202 self.firstdatatime = self.dataIn.utctime
210 self.firstdatatime = self.dataIn.utctime
203
211
204 if self.profIndex == nProfiles:
212 if self.profIndex == nProfiles or zeroPad:
205
213
206 self.__updateSpecFromVoltage()
214 self.__updateSpecFromVoltage()
207
215
@@ -682,6 +690,7 class getNoiseB(Operation):
682 #print(self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.nIncohInt)
690 #print(self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.nIncohInt)
683 noise = numpy.zeros( self.dataOut.nChannels)
691 noise = numpy.zeros( self.dataOut.nChannels)
684 norm = 1
692 norm = 1
693
685 for channel in range( self.dataOut.nChannels):
694 for channel in range( self.dataOut.nChannels):
686 if not hasattr(self.dataOut.nIncohInt,'__len__'):
695 if not hasattr(self.dataOut.nIncohInt,'__len__'):
687 norm = 1
696 norm = 1
@@ -691,10 +700,11 class getNoiseB(Operation):
691 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
700 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
692 daux = numpy.multiply(daux, norm)
701 daux = numpy.multiply(daux, norm)
693 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
702 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
694 #noise[channel] = self.getNoiseByMean(daux)/self.offset
703 # noise[channel] = self.getNoiseByMean(daux)/self.offset
695 #print(daux.shape, daux)
704 #print(daux.shape, daux)
696 noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
705 #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
697
706 sortdata = numpy.sort(daux, axis=None)
707 noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset
698 # data = numpy.mean(daux,axis=1)
708 # data = numpy.mean(daux,axis=1)
699 # sortdata = numpy.sort(data, axis=None)
709 # sortdata = numpy.sort(data, axis=None)
700 # noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset
710 # noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset
@@ -743,7 +753,7 class getNoiseB(Operation):
743 sortdata = numpy.sort(data, axis=None)
753 sortdata = numpy.sort(data, axis=None)
744
754
745 lenOfData = len(sortdata)
755 lenOfData = len(sortdata)
746 nums_min = lenOfData*0.05
756 nums_min = lenOfData*0.2
747
757
748 if nums_min <= 5:
758 if nums_min <= 5:
749
759
@@ -761,7 +771,7 class getNoiseB(Operation):
761 sumq += sortdata[j]**2
771 sumq += sortdata[j]**2
762 #sumq -= sump**2
772 #sumq -= sump**2
763 if j > nums_min:
773 if j > nums_min:
764 rtest = float(j)/(j-1) + 1.0/0.1
774 rtest = float(j)/(j-1) + 1.0/navg
765 #if ((sumq*j) > (sump**2)):
775 #if ((sumq*j) > (sump**2)):
766 if ((sumq*j) > (rtest*sump**2)):
776 if ((sumq*j) > (rtest*sump**2)):
767 j = j - 1
777 j = j - 1
@@ -1648,8 +1658,74 class IntegrationFaradaySpectra(Operation):
1648 #print(self.dataOut.flagNoData)
1658 #print(self.dataOut.flagNoData)
1649 return self.dataOut
1659 return self.dataOut
1650
1660
1661
1662
1651 class removeInterference(Operation):
1663 class removeInterference(Operation):
1652
1664
1665 def removeInterference3(self, min_hei = None, max_hei = None):
1666
1667 jspectra = self.dataOut.data_spc
1668 #jcspectra = self.dataOut.data_cspc
1669 jnoise = self.dataOut.getNoise()
1670 num_incoh = self.dataOut.max_nIncohInt
1671 #print(jspectra.shape)
1672 num_channel, num_prof, num_hei = jspectra.shape
1673 minHei = min_hei
1674 maxHei = max_hei
1675 ########################################################################
1676 if minHei == None or (minHei < self.dataOut.heightList[0]):
1677 minHei = self.dataOut.heightList[0]
1678
1679 if maxHei == None or (maxHei > self.dataOut.heightList[-1]):
1680 maxHei = self.dataOut.heightList[-1]
1681 minIndex = 0
1682 maxIndex = 0
1683 heights = self.dataOut.heightList
1684
1685 inda = numpy.where(heights >= minHei)
1686 indb = numpy.where(heights <= maxHei)
1687
1688 try:
1689 minIndex = inda[0][0]
1690 except:
1691 minIndex = 0
1692 try:
1693 maxIndex = indb[0][-1]
1694 except:
1695 maxIndex = len(heights)
1696
1697 if (minIndex < 0) or (minIndex > maxIndex):
1698 raise ValueError("some value in (%d,%d) is not valid" % (
1699 minIndex, maxIndex))
1700 if (maxIndex >= self.dataOut.nHeights):
1701 maxIndex = self.dataOut.nHeights - 1
1702
1703 ########################################################################
1704
1705
1706 #dataOut.max_nIncohInt * dataOut.nCohInt
1707 norm = self.dataOut.nIncohInt /self.dataOut.max_nIncohInt
1708 #print(norm.shape)
1709 # Subrutina de Remocion de la Interferencia
1710 for ich in range(num_channel):
1711 # Se ordena los espectros segun su potencia (menor a mayor)
1712 #power = jspectra[ich, mask_prof, :]
1713 interf = jspectra[ich, :, minIndex:maxIndex]
1714 #print(interf.shape)
1715 inttef = interf.mean(axis=1)
1716
1717 for hei in range(num_hei):
1718 temp = jspectra[ich,:, hei]
1719 temp -= inttef
1720 temp += jnoise[ich]*norm[ich,hei]
1721 jspectra[ich,:, hei] = temp
1722
1723 # Guardar Resultados
1724 self.dataOut.data_spc = jspectra
1725 #self.dataOut.data_cspc = jcspectra
1726
1727 return 1
1728
1653 def removeInterference2(self):
1729 def removeInterference2(self):
1654
1730
1655 cspc = self.dataOut.data_cspc
1731 cspc = self.dataOut.data_cspc
@@ -1687,7 +1763,7 class removeInterference(Operation):
1687
1763
1688 # hei_interf
1764 # hei_interf
1689 if hei_interf is None:
1765 if hei_interf is None:
1690 count_hei = int(num_hei / 2)
1766 count_hei = int(num_hei / 2) # a half of total ranges
1691 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1767 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1692 hei_interf = numpy.asarray(hei_interf)[0]
1768 hei_interf = numpy.asarray(hei_interf)[0]
1693 #print(hei_interf)
1769 #print(hei_interf)
@@ -1720,7 +1796,7 class removeInterference(Operation):
1720 power = power[:, hei_interf]
1796 power = power[:, hei_interf]
1721 power = power.sum(axis=0)
1797 power = power.sum(axis=0)
1722 psort = power.ravel().argsort()
1798 psort = power.ravel().argsort()
1723
1799 print(hei_interf[psort[list(range(offhei_interf, nhei_interf + offhei_interf))]])
1724 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1800 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1725 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1801 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1726 offhei_interf, nhei_interf + offhei_interf))]]]
1802 offhei_interf, nhei_interf + offhei_interf))]]]
@@ -1730,7 +1806,7 class removeInterference(Operation):
1730 tmp_noise = jnoise[ich]
1806 tmp_noise = jnoise[ich]
1731 junkspc_interf = junkspc_interf - tmp_noise
1807 junkspc_interf = junkspc_interf - tmp_noise
1732 #junkspc_interf[:,comp_mask_prof] = 0
1808 #junkspc_interf[:,comp_mask_prof] = 0
1733
1809 print(junkspc_interf.shape)
1734 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1810 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1735 jspc_interf = jspc_interf.transpose()
1811 jspc_interf = jspc_interf.transpose()
1736 # Calculando el espectro de interferencia promedio
1812 # Calculando el espectro de interferencia promedio
@@ -1743,7 +1819,6 class removeInterference(Operation):
1743
1819
1744 if (cnoiseid > 0):
1820 if (cnoiseid > 0):
1745 jspc_interf[noiseid] = 0
1821 jspc_interf[noiseid] = 0
1746
1747 # Expandiendo los perfiles a limpiar
1822 # Expandiendo los perfiles a limpiar
1748 if (cinterfid > 0):
1823 if (cinterfid > 0):
1749 new_interfid = (
1824 new_interfid = (
@@ -1757,16 +1832,14 class removeInterference(Operation):
1757
1832
1758 for ip in range(new_cinterfid):
1833 for ip in range(new_cinterfid):
1759 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1834 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1760 jspc_interf[new_interfid[ip]
1835 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1761 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1762
1836
1763 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1837 jspectra[ich, :, ind_hei] = jspectra[ich, :,ind_hei] - jspc_interf # Corregir indices
1764 ind_hei] - jspc_interf # Corregir indices
1765
1838
1766 # Removiendo la interferencia del punto de mayor interferencia
1839 # Removiendo la interferencia del punto de mayor interferencia
1767 ListAux = jspc_interf[mask_prof].tolist()
1840 ListAux = jspc_interf[mask_prof].tolist()
1768 maxid = ListAux.index(max(ListAux))
1841 maxid = ListAux.index(max(ListAux))
1769
1842 print(cinterfid)
1770 if cinterfid > 0:
1843 if cinterfid > 0:
1771 for ip in range(cinterfid * (interf == 2) - 1):
1844 for ip in range(cinterfid * (interf == 2) - 1):
1772 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1845 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
@@ -1783,16 +1856,15 class removeInterference(Operation):
1783
1856
1784 for id1 in range(4):
1857 for id1 in range(4):
1785 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1858 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1786
1787 xx_inv = numpy.linalg.inv(xx)
1859 xx_inv = numpy.linalg.inv(xx)
1788 xx = xx_inv[:, 0]
1860 xx = xx_inv[:, 0]
1789 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1861 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1790 yy = jspectra[ich, mask_prof[ind], :]
1862 yy = jspectra[ich, mask_prof[ind], :]
1791 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1863 jspectra[ich, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1792 yy.transpose(), xx)
1793
1864
1794 indAux = (jspectra[ich, :, :] < tmp_noise *
1865 indAux = (jspectra[ich, :, :] < tmp_noise *
1795 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1866 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1867 print(indAux)
1796 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1868 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1797 (1 - 1 / numpy.sqrt(num_incoh))
1869 (1 - 1 / numpy.sqrt(num_incoh))
1798
1870
@@ -1856,7 +1928,7 class removeInterference(Operation):
1856
1928
1857 return 1
1929 return 1
1858
1930
1859 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1931 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1, minHei=None, maxHei=None):
1860
1932
1861 self.dataOut = dataOut
1933 self.dataOut = dataOut
1862
1934
@@ -1864,7 +1936,8 class removeInterference(Operation):
1864 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1936 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1865 elif mode == 2:
1937 elif mode == 2:
1866 self.removeInterference2()
1938 self.removeInterference2()
1867
1939 elif mode == 3:
1940 self.removeInterference3(min_hei=minHei, max_hei=maxHei)
1868 return self.dataOut
1941 return self.dataOut
1869
1942
1870
1943
@@ -2051,7 +2124,7 class IncohInt(Operation):
2051 dataOut.data_outlier = self.nOutliers
2124 dataOut.data_outlier = self.nOutliers
2052 dataOut.utctime = avgdatatime
2125 dataOut.utctime = avgdatatime
2053 dataOut.flagNoData = False
2126 dataOut.flagNoData = False
2054 dataOut.max_nIncohInt += self.__profIndex
2127 dataOut.max_nIncohInt = self.__profIndex
2055 self.incohInt = 0
2128 self.incohInt = 0
2056 self.nOutliers = 0
2129 self.nOutliers = 0
2057 self.__profIndex = 0
2130 self.__profIndex = 0
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (588 lines changed) Show them Hide them
@@ -6,7 +6,7 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 from schainpy.utils import log
6 from schainpy.utils import log
7 from schainpy.model.io.utils import getHei_index
7 from schainpy.model.io.utils import getHei_index
8 from time import time
8 from time import time
9
9 import datetime
10 import numpy
10 import numpy
11 #import copy
11 #import copy
12 from schainpy.model.data import _noise
12 from schainpy.model.data import _noise
@@ -413,6 +413,42 class printAttribute(Operation):
413 if hasattr(dataOut, attr):
413 if hasattr(dataOut, attr):
414 log.log(getattr(dataOut, attr), attr)
414 log.log(getattr(dataOut, attr), attr)
415
415
416 class cleanHeightsInterf(Operation):
417
418 def run(self, dataOut, heightsList, repeats=0, step=0, factor=1, idate=None, startH=None, endH=None):
419
420 #print(dataOut.data.shape)
421
422 startTime = datetime.datetime.combine(idate,startH)
423 endTime = datetime.datetime.combine(idate,endH)
424 currentTime = datetime.datetime.fromtimestamp(dataOut.utctime)
425
426 if currentTime < startTime or currentTime > endTime:
427 return dataOut
428
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
439 for ch in range(dataOut.data.shape[0]):
440 i = 0
441 for hei in heights_indx:
442 #print(dataOut.data[ch,hei])
443 if dataOut.data.ndim < 3:
444 dataOut.data[ch,hei-1] = (dataOut.data[ch,hei-1])*wMask[i]
445 else:
446 dataOut.data[ch,:,hei-1] = (dataOut.data[ch,:,hei-1])*wMask[i]
447 #print("done")
448 i += 1
449
450 return dataOut
451
416
452
417 class interpolateHeights(Operation):
453 class interpolateHeights(Operation):
418
454
@@ -1785,7 +1821,7 class SSheightProfiles2(Operation):
1785
1821
1786 residue = (self.__nHeis - self.nsamples) % self.step
1822 residue = (self.__nHeis - self.nsamples) % self.step
1787 if residue != 0:
1823 if residue != 0:
1788 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
1824 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,self.__nProfiles - self.nsamples,residue))
1789
1825
1790 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1826 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1791 self.init_range = dataOut.heightList[0]
1827 self.init_range = dataOut.heightList[0]
@@ -1808,14 +1844,15 class SSheightProfiles2(Operation):
1808
1844
1809 if repeat is not None:
1845 if repeat is not None:
1810 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1846 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1811 if data.ndim == 2:
1847 if data.ndim < 3:
1812 data = data.reshape(1,1,self.__nHeis )
1848 data = data.reshape(self.__nChannels,1,self.__nHeis )
1813 #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape)
1849 #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape, code_block.shape)
1814 for i in range(int(self.new_nHeights)): #nuevas alturas
1850 for ch in range(self.__nChannels):
1815 if code is not None:
1851 for i in range(int(self.new_nHeights)): #nuevas alturas
1816 self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]*code_block
1852 if code is not None:
1817 else:
1853 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]*code_block
1818 self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1854 else:
1855 self.buffer[ch,i,:] = data[ch,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1819
1856
1820 for j in range(self.__nChannels): #en los cananles
1857 for j in range(self.__nChannels): #en los cananles
1821 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
1858 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
@@ -1824,7 +1861,7 class SSheightProfiles2(Operation):
1824
1861
1825
1862
1826 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1863 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1827
1864 # print("running")
1828 if dataOut.flagNoData == True:
1865 if dataOut.flagNoData == True:
1829 return dataOut
1866 return dataOut
1830 dataOut.flagNoData = True
1867 dataOut.flagNoData = True
@@ -1944,6 +1981,7 class removeProfileByFaradayHS(Operation):
1944
1981
1945 self.nChannels = dataOut.nChannels
1982 self.nChannels = dataOut.nChannels
1946 self.nHeights = dataOut.nHeights
1983 self.nHeights = dataOut.nHeights
1984 self.test_counter = 0
1947
1985
1948 def filterSatsProfiles(self):
1986 def filterSatsProfiles(self):
1949 data = self.__buffer_data
1987 data = self.__buffer_data
@@ -2020,6 +2058,137 class removeProfileByFaradayHS(Operation):
2020 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2058 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2021 return data
2059 return data
2022
2060
2061 def cleanSpikesFFT2D(self):
2062 incoh_int = 10
2063 norm_img = 75
2064 import matplotlib.pyplot as plt
2065 import datetime
2066 import cv2
2067 data = self.__buffer_data
2068 print("cleaning shape inpt: ",data.shape)
2069 self.__buffer_data = []
2070
2071
2072 channels , profiles, heights = data.shape
2073 len_split_prof = profiles / incoh_int
2074
2075
2076 for ch in range(channels):
2077 data_10 = numpy.split(data[ch, :, self.minHei_idx:], incoh_int, axis=0) # división de los perfiles
2078 print("splited data: ",len(data_10)," -> ", data_10[0].shape)
2079 int_img = None
2080 i_count = 0
2081 n_x, n_y = data_10[0].shape
2082 for s_data in data_10: #porciones de espectro
2083 spectrum = numpy.fft.fft2(s_data, axes=(0,1))
2084 z = numpy.abs(spectrum)
2085 mg = z[2:n_y,:] #omitir dc y adjunto
2086 dat = numpy.log10(mg.T)
2087 i_count += 1
2088 if i_count == 1:
2089 int_img = dat
2090 else:
2091 int_img += dat
2092 #print(i_count)
2093
2094 min, max = int_img.min(), int_img.max()
2095 int_img = ((int_img-min)*255/(max-min)).astype(numpy.uint8)
2096
2097 cv2.imshow('integrated image', int_img) #numpy.fft.fftshift(img))
2098 cv2.waitKey(0)
2099 #####################################################################
2100 kernel_h = numpy.zeros((3,3)) #
2101 kernel_h[0, :] = -2
2102 kernel_h[1, :] = 3
2103 kernel_h[2, :] = -2
2104
2105
2106 kernel_5h = numpy.zeros((5,5)) #
2107 kernel_5h[0, :] = -2
2108 kernel_5h[1, :] = -1
2109 kernel_5h[2, :] = 5
2110 kernel_5h[3, :] = -1
2111 kernel_5h[4, :] = -2
2112
2113 #####################################################################
2114 sharp_img = cv2.filter2D(src=int_img, ddepth=-1, kernel=kernel_5h)
2115 # cv2.imshow('sharp image h ', sharp_img)
2116 # cv2.waitKey(0)
2117 #####################################################################
2118 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (5,1)) #11
2119 #####################################################################
2120 detected_lines_h = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=1)
2121 # cv2.imshow('lines horizontal', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2122 # cv2.waitKey(0)
2123 #####################################################################
2124 ret, detected_lines_h = cv2.threshold(detected_lines_h, 200, 255, cv2.THRESH_BINARY)#
2125 cv2.imshow('binary img', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2126 cv2.waitKey(0)
2127 #####################################################################
2128 cnts_h, h0 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2129 #####################################################################
2130 h_line_index = []
2131 v_line_index = []
2132
2133 #cnts_h += cnts_h_s #combine large and small lines
2134
2135 # line indexes x1, x2, y
2136 for c in cnts_h:
2137 #print(c)
2138 if len(c) < 3: #contorno linea
2139 x1 = c[0][0][0]
2140 x2 = c[1][0][0]
2141 if x1 > 5 and x2 < (n_x-5) :
2142 start = incoh_int + (x1 * incoh_int)
2143 end = incoh_int + (x2 * incoh_int)
2144 h_line_index.append( [start, end, c[0][0][1]] )
2145
2146 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2147 else: #contorno poligono
2148 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2149 y = numpy.unique(pairs[:,1])
2150 x = numpy.unique(pairs[:,0])
2151 #print(x)
2152 for yk in y:
2153 x0 = x[0]
2154 if x0 < 8:
2155 x0 = 10
2156 #print(x[0], x[-1], yk)
2157 h_line_index.append( [x0, x[-1], yk])
2158 #print("x1, x2, y ->p ", x[0], x[-1], yk)
2159 ###################################################################
2160 #print("Cleaning")
2161 # # clean Spectrum
2162 spectrum = numpy.fft.fft2(data[ch,:,self.minHei_idx:], axes=(0,1))
2163 z = numpy.abs(spectrum)
2164 phase = numpy.angle(spectrum)
2165 print("Total Horizontal", len(h_line_index))
2166 if len(h_line_index) < 75 :
2167 for x1, x2, y in h_line_index:
2168 print(x1, x2, y)
2169 z[x1:x2,y] = 0
2170
2171
2172 spcCleaned = z * numpy.exp(1j*phase)
2173
2174 dat2 = numpy.log10(z[1:-1,:].T)
2175 min, max =dat2.min(), dat2.max()
2176 print(min, max)
2177 img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
2178 cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned))
2179 cv2.waitKey(0)
2180 cv2.destroyAllWindows()
2181
2182 data[ch,:,self.minHei_idx:] = numpy.fft.ifft2(spcCleaned, axes=(0,1))
2183
2184
2185 #print("cleanOutliersByBlock Done", data.shape)
2186 self.__buffer_data = data
2187 return data
2188
2189
2190
2191
2023 def cleanOutliersByBlock(self):
2192 def cleanOutliersByBlock(self):
2024 import matplotlib.pyplot as plt
2193 import matplotlib.pyplot as plt
2025 import datetime
2194 import datetime
@@ -2027,10 +2196,10 class removeProfileByFaradayHS(Operation):
2027 #print(self.__buffer_data[0].shape)
2196 #print(self.__buffer_data[0].shape)
2028 data = self.__buffer_data#.copy()
2197 data = self.__buffer_data#.copy()
2029 print("cleaning shape inpt: ",data.shape)
2198 print("cleaning shape inpt: ",data.shape)
2030
2031 self.__buffer_data = []
2199 self.__buffer_data = []
2032
2200
2033 spectrum = numpy.fft.fft2(data[:,:,self.minHei_idx:], axes=(0,2))
2201
2202 spectrum = numpy.fft.fft2(data[:,:,self.minHei_idx:], axes=(1,2))
2034 print("spc : ",spectrum.shape)
2203 print("spc : ",spectrum.shape)
2035 (nch,nsamples, nh) = spectrum.shape
2204 (nch,nsamples, nh) = spectrum.shape
2036 data2 = None
2205 data2 = None
@@ -2045,134 +2214,340 class removeProfileByFaradayHS(Operation):
2045 freqv = numpy.fft.fftfreq(nh, d=dt1)
2214 freqv = numpy.fft.fftfreq(nh, d=dt1)
2046 freqh = numpy.fft.fftfreq(self.n, d=dt2)
2215 freqh = numpy.fft.fftfreq(self.n, d=dt2)
2047
2216
2048 # freqv = numpy.fft.fftshift(numpy.fft.fftfreq(nh, d=dt1))
2217 z = numpy.abs(spectrum[ch,:,:])
2049 # freqh = numpy.fft.fftshift(numpy.fft.fftfreq(self.n, d=dt2))
2218 phase = numpy.angle(spectrum[ch,:,:])
2050 #
2219 z1 = z[0,:]
2220 #print("shape z: ", z.shape, nsamples)
2051
2221
2222 dat = numpy.log10(z[1:nsamples,:].T)
2052
2223
2224 pdat = numpy.log10(phase.T)
2225 #print("dat mean",dat.mean())
2053
2226
2227 mean, min, max = dat.mean(), dat.min(), dat.max()
2228 img = ((dat-min)*200/(max-min)).astype(numpy.uint8)
2054
2229
2230 # print(img.shape)
2231 cv2.imshow('image', img) #numpy.fft.fftshift(img))
2232 cv2.waitKey(0)
2055
2233
2056
2234
2057 z = numpy.abs(spectrum[ch,:,:])
2235 ''' #FUNCIONA LINEAS PEQUEÑAS
2058 phase = numpy.angle(spectrum[ch,:,:])
2236 kernel_5h = numpy.zeros((5,3)) #
2237 kernel_5h[0, :] = 2
2238 kernel_5h[1, :] = 1
2239 kernel_5h[2, :] = 0
2240 kernel_5h[3, :] = -1
2241 kernel_5h[4, :] = -2
2059
2242
2060 print("shape z: ", z.shape)
2243 sharp_imgh = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h)
2244 cv2.imshow('sharp image h',sharp_imgh)
2245 cv2.waitKey(0)
2246 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (20,1))
2061
2247
2248 detected_lines_h = cv2.morphologyEx(sharp_imgh, cv2.MORPH_OPEN, horizontal_kernel, iterations=1)
2249 #detected_lines_h = cv2.medianBlur(detected_lines_h, 3)
2250 #detected_lines_h = cv2.filter2D(src=img, ddepth=-1, kernel=kernel)
2251 cv2.imshow('lines h gray', detected_lines_h)
2252 cv2.waitKey(0)
2253 reth, detected_lines_h = cv2.threshold(detected_lines_h, 90, 255, cv2.THRESH_BINARY)
2254 cv2.imshow('lines h ', detected_lines_h)
2255 cv2.waitKey(0)
2256 '''
2062
2257
2063 # Find all peaks higher than the 98th percentile
2064 peaks = z < numpy.percentile(z, 99)
2065 #print(peaks)
2066 # Set those peak coefficients to zero
2067 #spectrum2 = spectrum * peaks.astype(int)
2068 #data2 = numpy.fft.ifft2(spectrum2,axes=(0,2))
2069
2258
2070 dat = numpy.log10(z.T)
2259 '''
2071 pdat = numpy.log10(phase.T)
2260 kernel_3h = numpy.zeros((3,10)) #10
2261 kernel_3h[0, :] = -1
2262 kernel_3h[1, :] = 2
2263 kernel_3h[2, :] = -1
2072
2264
2073
2265
2074 min, max = dat.min(), dat.max()
2266 kernel_h = numpy.zeros((3,20)) #20
2075 #img = ((dat-min)*255/(max-min)).astype(numpy.uint8)
2267 kernel_h[0, :] = -1
2076 img = ((dat-min)*200/(max-min)).astype(numpy.uint8)
2268 kernel_h[1, :] = 2
2077 #img = 255 - ((dat-min)*255/max).astype(numpy.uint8)
2269 kernel_h[2, :] = -1
2078
2270
2079 #print(img.shape)
2271 kernel_v = numpy.zeros((30,3)) #30
2080 # cv2.imshow('image', numpy.fft.fftshift(img))
2272 kernel_v[:, 0] = -1
2081 # cv2.waitKey(0)
2273 kernel_v[:, 1] = 2
2274 kernel_v[:, 2] = -1
2275
2276 kernel_4h = numpy.zeros((4,20)) #
2277 kernel_4h[0, :] = 1
2278 kernel_4h[1, :] = 0
2279 kernel_4h[2, :] = 0
2280 kernel_4h[3, :] = -1
2082
2281
2083 spikes = numpy.where(img > 230, True, False)
2282 kernel_5h = numpy.zeros((5,30)) #
2084 #print("spikes: ", spikes.sum())
2283 kernel_5h[0, :] = 2
2284 kernel_5h[1, :] = 1
2285 kernel_5h[2, :] = 0
2286 kernel_5h[3, :] = -1
2287 kernel_5h[4, :] = -2
2085
2288
2086 #img = cv2.medianBlur(img,3)
2087
2289
2088 #cv2.imshow('image', cv2.resize(img.astype('float32'), (600, 600)))
2290 sharp_img0 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_3h)
2089 kernel3 = numpy.zeros((3,15))
2291 # cv2.imshow('sharp image small h',sharp_img0) # numpy.fft.fftshift(sharp_img1))
2090 kernel3[0, :] = -1
2292 # cv2.waitKey(0)
2091 kernel3[1, :] = 2
2293
2092 kernel3[2, :] = -1
2294 sharp_img1 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_h)
2093 sharp_img = cv2.filter2D(src=img, ddepth=-1, kernel=kernel3)
2295 # cv2.imshow('sharp image h',sharp_img1) # numpy.fft.fftshift(sharp_img1))
2296 # cv2.waitKey(0)
2094
2297
2298 sharp_img2 = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_v)
2299 # cv2.imshow('sharp image v', sharp_img2) #numpy.fft.fftshift(sharp_img2))
2300 # cv2.waitKey(0)
2095
2301
2096 # cv2.imshow('sharp image', numpy.fft.fftshift(sharp_img))
2302 sharp_imgw = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_4h)
2303 # cv2.imshow('sharp image h wide', sharp_imgw) #numpy.fft.fftshift(sharp_img2))
2097 # cv2.waitKey(0)
2304 # cv2.waitKey(0)
2098
2305
2099 ret, thresh = cv2.threshold(sharp_img, 127, 255, cv2.THRESH_BINARY_INV + cv2.THRESH_OTSU)
2306 sharp_imgwl = cv2.filter2D(src=img, ddepth=-1, kernel=kernel_5h, borderType = cv2.BORDER_ISOLATED)
2100 # #thresh = cv2.adaptiveThreshold(sharp_img,255,cv2.ADAPTIVE_THRESH_MEAN_C,cv2.THRESH_BINARY_INV,21,0)
2307 cv2.imshow('sharp image h long wide', sharp_imgwl) #numpy.fft.fftshift(sharp_img2))
2101 # cv2.imshow('binary img', numpy.fft.fftshift(thresh))
2308 cv2.waitKey(0)
2309
2310 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_h.jpg', sharp_img1)
2311 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/sharp_v.jpg', sharp_img2)
2312 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/input_img.jpg', img)
2313
2314 ########################small horizontal
2315 horizontal_kernel_s = cv2.getStructuringElement(cv2.MORPH_RECT, (11,1)) #11
2316 ######################## horizontal
2317 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) #30
2318 ######################## vertical
2319 vertical_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (1,50)) #50
2320 ######################## horizontal wide
2321 horizontal_kernel_w = cv2.getStructuringElement(cv2.MORPH_RECT, (30,1)) # 30
2322
2323 horizontal_kernel_expand = cv2.getStructuringElement(cv2.MORPH_RECT, (3,3)) #
2324
2325 horizontal_kernel_wl = cv2.getStructuringElement(cv2.MORPH_RECT, (50,1)) #
2326
2327 detected_lines_h_s = cv2.morphologyEx(sharp_img0, cv2.MORPH_OPEN, horizontal_kernel_s, iterations=7) #7
2328 detected_lines_h = cv2.morphologyEx(sharp_img1, cv2.MORPH_OPEN, horizontal_kernel, iterations=7) #7
2329 detected_lines_v = cv2.morphologyEx(sharp_img2, cv2.MORPH_OPEN, vertical_kernel, iterations=7) #7
2330 detected_lines_h_w = cv2.morphologyEx(sharp_imgw, cv2.MORPH_OPEN, horizontal_kernel_w, iterations=5) #5
2331
2332 detected_lines_h_wl = cv2.morphologyEx(sharp_imgwl, cv2.MORPH_OPEN, horizontal_kernel_wl, iterations=5) #
2333 detected_lines_h_wl = cv2.filter2D(src=detected_lines_h_wl, ddepth=-1, kernel=horizontal_kernel_expand)
2334
2335 # cv2.imshow('lines h small gray', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h))
2336 # cv2.waitKey(0)
2337 # cv2.imshow('lines h gray', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2102 # cv2.waitKey(0)
2338 # cv2.waitKey(0)
2103 # #Remove horizontal
2339 # cv2.imshow('lines v gray', detected_lines_v) #numpy.fft.fftshift(detected_lines_h))
2104 horizontal_kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (10,1))
2340 # cv2.waitKey(0)
2105 detected_lines = cv2.morphologyEx(sharp_img, cv2.MORPH_OPEN, horizontal_kernel, iterations=8)
2341 # cv2.imshow('lines h wide gray', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h))
2106 cnts, h = cv2.findContours(detected_lines, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2342 # cv2.waitKey(0)
2343 cv2.imshow('lines h long wide gray', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_h))
2344 cv2.waitKey(0)
2345
2346 reth_s, detected_lines_h_s = cv2.threshold(detected_lines_h_s, 85, 255, cv2.THRESH_BINARY)# 85
2347 reth, detected_lines_h = cv2.threshold(detected_lines_h, 30, 255, cv2.THRESH_BINARY) #30
2348 retv, detected_lines_v = cv2.threshold(detected_lines_v, 30, 255, cv2.THRESH_BINARY) #30
2349 reth_w, detected_lines_h_w = cv2.threshold(detected_lines_h_w, 35, 255, cv2.THRESH_BINARY)#
2350 reth_wl, detected_lines_h_wl = cv2.threshold(detected_lines_h_wl, 200, 255, cv2.THRESH_BINARY)#
2351
2352 cnts_h_s, h0 = cv2.findContours(detected_lines_h_s, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2353 cnts_h, h1 = cv2.findContours(detected_lines_h, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2354 cnts_v, h2 = cv2.findContours(detected_lines_v, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2355 cnts_h_w, h3 = cv2.findContours(detected_lines_h_w, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2356 cnts_h_wl, h4 = cv2.findContours(detected_lines_h_wl, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
2357 #print("horizontal ", cnts_h)
2358 #print("vertical ", cnts_v)
2359 # cnts, h = cv2.findContours(detected_lines_h, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
2107 # print(cnts)
2360 # print(cnts)
2108 # cv2.imshow('lines ', numpy.fft.fftshift(detected_lines))
2361 # cv2.imshow('lines h wide', detected_lines_h_w) #numpy.fft.fftshift(detected_lines_h))
2362 # cv2.waitKey(0)
2363 cv2.imshow('lines h wide long ', detected_lines_h_wl) #numpy.fft.fftshift(detected_lines_v))
2364 cv2.waitKey(0)
2365 # cv2.imshow('lines h small', detected_lines_h_s) #numpy.fft.fftshift(detected_lines_h))
2366 # cv2.waitKey(0)
2367 # cv2.imshow('lines h ', detected_lines_h) #numpy.fft.fftshift(detected_lines_h))
2109 # cv2.waitKey(0)
2368 # cv2.waitKey(0)
2369 # cv2.imshow('lines v ', detected_lines_v) #numpy.fft.fftshift(detected_lines_v))
2370 # cv2.waitKey(0)
2371
2372 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_h.jpg', detected_lines_h)
2373 # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/lines_v.jpg', detected_lines_v)
2374
2110 #cnts = cnts[0] if len(cnts) == 2 else cnts[1]
2375 #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 ])
2376 #y_line_index = numpy.asarray([ [c[0][0][0],c[1][0][0], c[0][0][1]] for c in cnts_v ])
2112
2377 h_line_index = []
2113 # print(y_line_index)
2378 v_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
2379
2134 readed_lines.append(y)
2380 cnts_h += cnts_h_s #combine large and small lines
2381
2382 # line indexes x1, x2, y
2383 for c in cnts_h:
2384 #print(c)
2385 if len(c) < 3: #contorno linea
2386 x1 = c[0][0][0]
2387 if x1 < 8:
2388 x1 = 10
2389 h_line_index.append( [x1,c[1][0][0], c[0][0][1]] )
2390 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2391 else: #contorno poligono
2392 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2393 y = numpy.unique(pairs[:,1])
2394 x = numpy.unique(pairs[:,0])
2395 #print(x)
2396 for yk in y:
2397 x0 = x[0]
2398 if x0 < 8:
2399 x0 = 10
2400 #print(x[0], x[-1], yk)
2401 h_line_index.append( [x0, x[-1], yk])
2402 #print("x1, x2, y ->p ", x[0], x[-1], yk)
2403 for c in cnts_h_w:
2404 #print(c)
2405 if len(c) < 3: #contorno linea
2406 x1 = c[0][0][0]
2407 if x1 < 8:
2408 x1 = 10
2409 y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro
2410 h_line_index.append( [x1,c[1][0][0],y] )
2411 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2412 else: #contorno poligono
2413 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2414 y = numpy.unique(pairs[:,1])
2415 x = numpy.unique(pairs[:,0])
2416 #print(x)
2417 for yk in y:
2418
2419 x0 = x[0]
2420 if x0 < 8:
2421 x0 = 10
2422 h_line_index.append( [x0, x[-1], yk-2])
2423
2424 for c in cnts_h_wl: # # revisar
2425 #print(c)
2426 if len(c) < 3: #contorno linea
2427 x1 = c[0][0][0]
2428 if x1 < 8:
2429 x1 = 10
2430 y = c[0][0][1] - 2 # se incrementa 2 líneas x el filtro
2431 h_line_index.append( [x1,c[1][0][0],y] )
2432 #print("x1, x2, y", c[0][0][0],c[1][0][0], c[0][0][1])
2433 else: #contorno poligono
2434 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2435 y = numpy.unique(pairs[:,1])
2436 x = numpy.unique(pairs[:,0])
2437 for yk in range(y[-1]-y[0]):
2438 y_k = yk +y[0]
2439
2440 x0 = x[0]
2441 if x0 < 8:
2442 x0 = 10
2443 h_line_index.append( [x0, x[-1], y_k-2])
2444
2445 print([[c[0][0][1],c[1][0][1], c[0][0][0] ] for c in cnts_v])
2446 # line indexes y1, y2, x
2447 for c in cnts_v:
2448 if len(c) < 3: #contorno linea
2449 v_line_index.append( [c[0][0][1],c[1][0][1], c[0][0][0] ] )
2450 else: #contorno poligono
2451 pairs = numpy.asarray([c[n][0] for n in range(len(c))])
2452 #print(pairs)
2453 y = numpy.unique(pairs[:,1])
2454 x = numpy.unique(pairs[:,0])
2455
2456 for xk in x:
2457 #print(x[0], x[-1], yk)
2458 v_line_index.append( [y[0],y[-1], xk])
2459
2460 ###################################################################
2461 # # clean Horizontal
2462 print("Total Horizontal", len(h_line_index))
2463 if len(h_line_index) < 75 :
2464 for x1, x2, y in h_line_index:
2465 #print("cleaning ",x1, x2, y)
2466 len_line = x2 - x1
2467 if y > 10 and y < (nh -10):
2468 # if y != (nh-1):
2469 # list = [ ((z[n, y-1] + z[n,y+1])/2) for n in range(len_line)]
2470 # else:
2471 # list = [ ((z[n, y-1] + z[n,0])/2) for n in range(len_line)]
2472 #
2473 # z[x1:x2,y] = numpy.asarray(list)
2474 z[x1-5:x2+5,y:y+1] = 0
2475
2476 # clean vertical
2477 for y1, y2, x in v_line_index:
2478 len_line = y2 - y1
2479 #print(x)
2480 if x > 0 and x < (nsamples-2):
2481 # if x != (nsamples-1):
2482 # list = [ ((z[x-2, n] + z[x+2,n])/2) for n in range(len_line)]
2483 # else:
2484 # list = [ ((z[x-2, n] + z[1,n])/2) for n in range(len_line)]
2485 #
2486 # #z[x-1:x+1,y1:y2] = numpy.asarray(list)
2487 #
2488 z[x+1,y1:y2] = 0
2489
2490 '''
2491 #z[: ,[215, 217, 221, 223, 225, 340, 342, 346, 348, 350, 465, 467, 471, 473, 475]]=0
2492 z[1: ,[112, 114, 118, 120, 122, 237, 239, 245, 247, 249, 362, 364, 368, 370, 372]]=0
2493 # z[: ,217]=0
2494 # z[: ,221]=0
2495 # z[: ,223]=0
2496 # z[: ,225]=0
2135
2497
2136 dat2 = numpy.log10(z.T)
2498 dat2 = numpy.log10(z.T)
2137 img = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
2499 #print(dat2)
2500 max = dat2.max()
2501 #print(" min, max ", max, min)
2502 img2 = ((dat2-min)*255/(max-min)).astype(numpy.uint8)
2503 #img_cleaned = img2.copy()
2504 #cv2.drawContours(img2, cnts_h, -1, (255,255,255), 1)
2505 #cv2.drawContours(img2, cnts_v, -1, (255,255,255), 1)
2138
2506
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
2507
2144 spcCleaned = z * numpy.exp(1j*phase)
2508 spcCleaned = z * numpy.exp(1j*phase)
2145 #print(spcCleaned)
2509 #print(spcCleaned)
2146
2510
2147
2511
2148 cv2.imshow('image2', numpy.fft.fftshift(img))
2512 # cv2.imshow('image contours', img2) #numpy.fft.fftshift(img))
2513 # cv2.waitKey(0)
2514
2515 cv2.imshow('cleaned', img2) #numpy.fft.fftshift(img_cleaned))
2149 cv2.waitKey(0)
2516 cv2.waitKey(0)
2517 # # cv2.imwrite('/home/soporte/Data/AMISR14/ISR/spc/spc/cleaned_{}.jpg'.format(self.test_counter), img2)
2150 cv2.destroyAllWindows()
2518 cv2.destroyAllWindows()
2519 # self.test_counter += 1
2151
2520
2152 m = numpy.mean(dat)
2153 o = numpy.std(dat)
2154
2521
2155 fig, ax = plt.subplots(1,2,figsize=(12, 6))
2522 #print("DC difference " ,z1 - z[0,:])
2156 #X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
2157 X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv))
2158
2523
2159 colormap = 'jet'
2524 # m = numpy.mean(dat)
2160 #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o))
2525 # o = numpy.std(dat)
2161 c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 4.2, vmax = 5.0)
2526 # print("mean ", m, " std ", o)
2162 fig.colorbar(c, ax=ax[0])
2527 # fig, ax = plt.subplots(1,2,figsize=(12, 6))
2163
2528 # #X, Y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
2164
2529 # X, Y = numpy.meshgrid(numpy.fft.fftshift(freqh),numpy.fft.fftshift(freqv))
2165 #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
2530 #
2166 #date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
2531 # colormap = 'jet'
2167 #strftime('%Y-%m-%d %H:%M:%S')
2532 # #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o))
2168 #ax[0].set_title('Spectrum magnitude '+date_time)
2533 # #c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = 6.5, vmax = 6.8)
2169 #fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
2534 # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o))
2170 #print("aqui estoy2",dat2[:,:,0].shape)
2535 # fig.colorbar(c, ax=ax[0])
2171 c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = 4.2, vmax = 5.0)
2536 #
2172 #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1)
2537 #
2173 #print("aqui estoy3")
2538 # #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
2174 fig.colorbar(c, ax=ax[1])
2539 # #date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
2175 plt.show()
2540 # #strftime('%Y-%m-%d %H:%M:%S')
2541 # #ax[0].set_title('Spectrum magnitude '+date_time)
2542 # #fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
2543 # #print("aqui estoy2",dat2[:,:,0].shape)
2544 # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap, vmin = 4.2, vmax = 5.0)
2545 # c = ax[0].pcolormesh(X, Y, numpy.fft.fftshift(dat2), cmap =colormap, vmin = (m-2*o), vmax = (m+1.5*o))
2546 # #c = ax[1].pcolormesh(X, Y, numpy.fft.fftshift(pdat), cmap =colormap ) #, vmin = 0.0, vmax = 0.5)
2547 # #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1)
2548 # #print("aqui estoy3")
2549 # fig.colorbar(c, ax=ax[1])
2550 # plt.show()
2176
2551
2177 spectrum[ch,:,:] = spcCleaned
2552 spectrum[ch,:,:] = spcCleaned
2178
2553
@@ -2180,10 +2555,10 class removeProfileByFaradayHS(Operation):
2180
2555
2181
2556
2182
2557
2183 data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(0,2))
2558 data[:,:,self.minHei_idx:] = numpy.fft.ifft2(spectrum, axes=(1,2))
2184
2559
2185 #print("cleanOutliersByBlock Done", data.shape)
2560 #print("cleanOutliersByBlock Done", data.shape)
2186
2561 self.__buffer_data = data
2187 return data
2562 return data
2188
2563
2189
2564
@@ -2212,6 +2587,7 class removeProfileByFaradayHS(Operation):
2212 if self.__profIndex == self.n:
2587 if self.__profIndex == self.n:
2213 #print("apnd : ",data)
2588 #print("apnd : ",data)
2214 #dataBlock = self.cleanOutliersByBlock()
2589 #dataBlock = self.cleanOutliersByBlock()
2590 #dataBlock = self.cleanSpikesFFT2D()
2215 dataBlock = self.filterSatsProfiles()
2591 dataBlock = self.filterSatsProfiles()
2216 self.__dataReady = True
2592 self.__dataReady = True
2217
2593
@@ -2288,7 +2664,7 class removeProfileByFaradayHS(Operation):
2288
2664
2289 self.init_prof = self.end_prof
2665 self.init_prof = self.end_prof
2290 self.end_prof += self.lenProfileOut
2666 self.end_prof += self.lenProfileOut
2291 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2667 # #print("data release shape: ",dataOut.data.shape, self.end_prof, dataOut.flagNoData)
2292 self.n_prof_released += 1
2668 self.n_prof_released += 1
2293
2669
2294 if self.end_prof >= (self.n +self.lenProfileOut):
2670 if self.end_prof >= (self.n +self.lenProfileOut):
General Comments 0
You need to be logged in to leave comments. Login now