##// END OF EJS Templates
schain
RSS
Eliminación de perfiles, retención de bloques, salida en pequeñps bloques o perfiles, rti de outliers
joabAM -
r1528:9eee9f011a2c
parent child
Show More
Side by Side Unified
Context file:
r1528:9eee9f011a2c Loading diff...:
Show file before | Show file after | Hide comments
@@ -66,7 +66,7 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
66 }
66 }
67 double *sortdata = (double*)PyArray_DATA(data_array);
67 double *sortdata = (double*)PyArray_DATA(data_array);
68 int lenOfData = (int)PyArray_SIZE(data_array) ;
68 int lenOfData = (int)PyArray_SIZE(data_array) ;
69 double nums_min = lenOfData*0.75;
69 double nums_min = lenOfData*0.85;
70 if (nums_min <= 5) nums_min = 5;
70 if (nums_min <= 5) nums_min = 5;
71 double sump = 0;
71 double sump = 0;
72 double sumq = 0;
72 double sumq = 0;
Show file before | Show file after | Hide comments
@@ -661,12 +661,10 class Project(Process):
661 elif ok == 'no_Read' and (not flag_no_read):
661 elif ok == 'no_Read' and (not flag_no_read):
662 nProc_noRead = n_proc - 1
662 nProc_noRead = n_proc - 1
663 flag_no_read = True
663 flag_no_read = True
664 print("No read")
665 continue
664 continue
666 elif ok == 'new_Read':
665 elif ok == 'new_Read':
667 nProc_noRead = 0
666 nProc_noRead = 0
668 flag_no_read = False
667 flag_no_read = False
669 print("read again")
670 continue
668 continue
671 elif not ok:
669 elif not ok:
672 #print("not ok",ok)
670 #print("not ok",ok)
Show file before | Show file after | Hide comments
@@ -78,7 +78,7 def hildebrand_sekhon(data, navg):
78 sortdata = numpy.sort(data, axis=None)
78 sortdata = numpy.sort(data, axis=None)
79 '''
79 '''
80 lenOfData = len(sortdata)
80 lenOfData = len(sortdata)
81 nums_min = lenOfData*0.2
81 nums_min = lenOfData*0.5
82
82
83 if nums_min <= 5:
83 if nums_min <= 5:
84
84
@@ -106,6 +106,7 def hildebrand_sekhon(data, navg):
106 j += 1
106 j += 1
107
107
108 lnoise = sump / j
108 lnoise = sump / j
109 return lnoise
109 '''
110 '''
110 return _noise.hildebrand_sekhon(sortdata, navg)
111 return _noise.hildebrand_sekhon(sortdata, navg)
111
112
@@ -380,7 +381,7 class Voltage(JROData):
380 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
381 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
381 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
382 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
382
383
383 def getNoisebyHildebrand(self, channel=None):
384 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
384 """
385 """
385 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
386 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
386
387
@@ -389,10 +390,10 class Voltage(JROData):
389 """
390 """
390
391
391 if channel != None:
392 if channel != None:
392 data = self.data[channel]
393 data = self.data[channel,ymin_index:ymax_index]
393 nChannels = 1
394 nChannels = 1
394 else:
395 else:
395 data = self.data
396 data = self.data[:,ymin_index:ymax_index]
396 nChannels = self.nChannels
397 nChannels = self.nChannels
397
398
398 noise = numpy.zeros(nChannels)
399 noise = numpy.zeros(nChannels)
@@ -407,10 +408,10 class Voltage(JROData):
407
408
408 return noise
409 return noise
409
410
410 def getNoise(self, type=1, channel=None):
411 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
411
412
412 if type == 1:
413 if type == 1:
413 noise = self.getNoisebyHildebrand(channel)
414 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
414
415
415 return noise
416 return noise
416
417
@@ -437,6 +438,8 class Voltage(JROData):
437
438
438 class Spectra(JROData):
439 class Spectra(JROData):
439
440
441 data_outlier = 0
442
440 def __init__(self):
443 def __init__(self):
441 '''
444 '''
442 Constructor
445 Constructor
@@ -475,6 +478,9 class Spectra(JROData):
475 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
478 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
476
479
477
480
481 self.max_nIncohInt = 1
482
483
478 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
484 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
479 """
485 """
480 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
486 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
@@ -482,12 +488,28 class Spectra(JROData):
482 Return:
488 Return:
483 noiselevel
489 noiselevel
484 """
490 """
485
491 # if hasattr(self.nIncohInt, "__len__"): #nIncohInt is a matrix
492 #
493 # heis = self.data_spc.shape[2]
494 #
495 # noise = numpy.zeros((self.nChannels, heis))
496 # for hei in range(heis):
497 # for channel in range(self.nChannels):
498 # daux = self.data_spc[channel, xmin_index:xmax_index, hei]
499 #
500 # noise[channel,hei] = hildebrand_sekhon(daux, self.nIncohInt[channel,hei])
501 #
502 # else:
503 # noise = numpy.zeros(self.nChannels)
504 # for channel in range(self.nChannels):
505 # daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
506 #
507 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
486 noise = numpy.zeros(self.nChannels)
508 noise = numpy.zeros(self.nChannels)
487 for channel in range(self.nChannels):
509 for channel in range(self.nChannels):
488 daux = self.data_spc[channel,
510 daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
489 xmin_index:xmax_index, ymin_index:ymax_index]
511
490 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
512 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt)
491
513
492 return noise
514 return noise
493
515
@@ -552,6 +574,7 class Spectra(JROData):
552 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
574 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
553 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
575 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
554
576
577
555 return normFactor
578 return normFactor
556
579
557 @property
580 @property
@@ -582,7 +605,7 class Spectra(JROData):
582 def getPower(self):
605 def getPower(self):
583
606
584 factor = self.normFactor
607 factor = self.normFactor
585 z = self.data_spc / factor
608 z = numpy.divide(self.data_spc,factor)
586 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
609 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
587 avg = numpy.average(z, axis=1)
610 avg = numpy.average(z, axis=1)
588
611
Show file before | Show file after | Hide comments
@@ -58,7 +58,7 class SpectraPlot(Plot):
58 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
58 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
59 data['spc'] = spc
59 data['spc'] = spc
60 data['rti'] = dataOut.getPower()
60 data['rti'] = dataOut.getPower()
61 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
61 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
62 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
62 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
63 if self.CODE == 'spc_moments':
63 if self.CODE == 'spc_moments':
64 data['moments'] = dataOut.moments
64 data['moments'] = dataOut.moments
@@ -86,9 +86,10 class SpectraPlot(Plot):
86
86
87 data = self.data[-1]
87 data = self.data[-1]
88 z = data['spc']
88 z = data['spc']
89 print(z.shape, x.shape, y.shape)
89 #print(z.shape, x.shape, y.shape)
90 for n, ax in enumerate(self.axes):
90 for n, ax in enumerate(self.axes):
91 noise = data['noise'][n]
91 #noise = data['noise'][n]
92 noise=62
92 if self.CODE == 'spc_moments':
93 if self.CODE == 'spc_moments':
93 mean = data['moments'][n, 1]
94 mean = data['moments'][n, 1]
94 if ax.firsttime:
95 if ax.firsttime:
@@ -105,15 +106,15 class SpectraPlot(Plot):
105 if self.showprofile:
106 if self.showprofile:
106 ax.plt_profile = self.pf_axes[n].plot(
107 ax.plt_profile = self.pf_axes[n].plot(
107 data['rti'][n], y)[0]
108 data['rti'][n], y)[0]
108 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
109 # ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
109 color="k", linestyle="dashed", lw=1)[0]
110 # color="k", linestyle="dashed", lw=1)[0]
110 if self.CODE == 'spc_moments':
111 if self.CODE == 'spc_moments':
111 ax.plt_mean = ax.plot(mean, y, color='k')[0]
112 ax.plt_mean = ax.plot(mean, y, color='k')[0]
112 else:
113 else:
113 ax.plt.set_array(z[n].T.ravel())
114 ax.plt.set_array(z[n].T.ravel())
114 if self.showprofile:
115 if self.showprofile:
115 ax.plt_profile.set_data(data['rti'][n], y)
116 ax.plt_profile.set_data(data['rti'][n], y)
116 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
117 #ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
117 if self.CODE == 'spc_moments':
118 if self.CODE == 'spc_moments':
118 ax.plt_mean.set_data(mean, y)
119 ax.plt_mean.set_data(mean, y)
119 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
120 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
@@ -274,7 +275,7 class RTIPlot(Plot):
274
275
275 self.x = self.data.times
276 self.x = self.data.times
276 self.y = self.data.yrange
277 self.y = self.data.yrange
277
278 #print(" x, y: ",self.x, self.y)
278 self.z = self.data[self.CODE]
279 self.z = self.data[self.CODE]
279 self.z = numpy.array(self.z, dtype=float)
280 self.z = numpy.array(self.z, dtype=float)
280 self.z = numpy.ma.masked_invalid(self.z)
281 self.z = numpy.ma.masked_invalid(self.z)
@@ -816,6 +817,7 class NoiselessSpectraPlot(Plot):
816 plot_type = 'pcolor'
817 plot_type = 'pcolor'
817 buffering = False
818 buffering = False
818 channelList = []
819 channelList = []
820 last_noise = None
819
821
820 def setup(self):
822 def setup(self):
821
823
@@ -842,18 +844,22 class NoiselessSpectraPlot(Plot):
842 self.update_list(dataOut)
844 self.update_list(dataOut)
843 data = {}
845 data = {}
844 meta = {}
846 meta = {}
845 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
846 (nch, nff, nh) = dataOut.data_spc.shape
847 n1 = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
848 noise = numpy.repeat(n1,nff, axis=1).reshape((nch,nff,nh))
849 #print(noise.shape, "noise", noise)
850
847
851 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - noise
848 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
849 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
852
850
853 data['spc'] = spc
851 if self.last_noise == None:
854 data['rti'] = dataOut.getPower() - n1
852 self.last_noise = n0
853 else:
854 n0 = (n0*0.2 + self.last_noise*0.8)
855 self.last_noise = n0
856
857 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
855
858
856 data['noise'] = n0
859 data['spc'] = spc - n0
860 data['rti'] = dataOut.getPower() - n0
861
862 #data['noise'] = noise
857 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
863 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
858
864
859 return data, meta
865 return data, meta
@@ -877,7 +883,7 class NoiselessSpectraPlot(Plot):
877 z = data['spc']
883 z = data['spc']
878
884
879 for n, ax in enumerate(self.axes):
885 for n, ax in enumerate(self.axes):
880 noise = data['noise'][n]
886 #noise = data['noise'][n]
881
887
882 if ax.firsttime:
888 if ax.firsttime:
883 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
889 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
@@ -893,16 +899,15 class NoiselessSpectraPlot(Plot):
893 if self.showprofile:
899 if self.showprofile:
894 ax.plt_profile = self.pf_axes[n].plot(
900 ax.plt_profile = self.pf_axes[n].plot(
895 data['rti'][n], y)[0]
901 data['rti'][n], y)[0]
896 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
902
897 color="k", linestyle="dashed", lw=1)[0]
898
903
899 else:
904 else:
900 ax.plt.set_array(z[n].T.ravel())
905 ax.plt.set_array(z[n].T.ravel())
901 if self.showprofile:
906 if self.showprofile:
902 ax.plt_profile.set_data(data['rti'][n], y)
907 ax.plt_profile.set_data(data['rti'][n], y)
903 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
904
908
905 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
909
910 self.titles.append('CH {}'.format(self.channelList[n]))
906
911
907
912
908 class NoiselessRTIPlot(Plot):
913 class NoiselessRTIPlot(Plot):
@@ -917,6 +922,7 class NoiselessRTIPlot(Plot):
917 channelList = []
922 channelList = []
918 elevationList = []
923 elevationList = []
919 azimuthList = []
924 azimuthList = []
925 last_noise = None
920
926
921 def setup(self):
927 def setup(self):
922 self.xaxis = 'time'
928 self.xaxis = 'time'
@@ -945,19 +951,21 class NoiselessRTIPlot(Plot):
945 data = {}
951 data = {}
946 meta = {}
952 meta = {}
947
953
948 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
954
949 (nch, nff, nh) = dataOut.data_spc.shape
955 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
950 #print(nch, nff, nh)
956 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
951 if nch != 1:
957 #print("noise: ",n0, dataOut.normFactor, norm, dataOut.nIncohInt, dataOut.max_nIncohInt)
952 aux = []
958 if self.last_noise == None:
953 for c in self.channelList:
959 self.last_noise = n0
954 aux.append(n0[c])
960 else:
955 n0 = numpy.asarray(aux)
961 n0 = (n0*0.2 + self.last_noise*0.8)
956 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
962 self.last_noise = n0
957 #print(dataOut.elevationList, dataOut.azimuthList)
963
958 #print(dataOut.channelList)
964
959 data['noiseless_rti'] = dataOut.getPower() - noise
965 data['noiseless_rti'] = dataOut.getPower() - n0
960 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
966
967 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
968 #print(noise)
961 return data, meta
969 return data, meta
962
970
963 def plot(self):
971 def plot(self):
@@ -986,6 +994,7 class NoiselessRTIPlot(Plot):
986 else:
994 else:
987 x, y, z = self.fill_gaps(*self.decimate())
995 x, y, z = self.fill_gaps(*self.decimate())
988 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
996 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
997 #print("plot shapes ", z.shape, x.shape, y.shape)
989 for n, ax in enumerate(self.axes):
998 for n, ax in enumerate(self.axes):
990 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
999 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
991 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1000 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
@@ -999,9 +1008,6 class NoiselessRTIPlot(Plot):
999 if self.showprofile:
1008 if self.showprofile:
1000 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1009 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1001
1010
1002 if "noise" in self.data:
1003 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
1004 color="k", linestyle="dashed", lw=1)[0]
1005 else:
1011 else:
1006 ax.collections.remove(ax.collections[0])
1012 ax.collections.remove(ax.collections[0])
1007 ax.plt = ax.pcolormesh(x, y, z[n].T,
1013 ax.plt = ax.pcolormesh(x, y, z[n].T,
@@ -1011,6 +1017,164 class NoiselessRTIPlot(Plot):
1011 )
1017 )
1012 if self.showprofile:
1018 if self.showprofile:
1013 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1019 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1014 if "noise" in self.data:
1020 # if "noise" in self.data:
1015 ax.plot_noise.set_data(numpy.repeat(
1021 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1016 data['noise'][n], len(self.y)), self.y)
1022 # ax.plot_noise.set_data(data['noise'][n], self.y)
1023
1024
1025 class OutliersRTIPlot(Plot):
1026 '''
1027 Plot for data_xxxx object
1028 '''
1029
1030 CODE = 'outlier'
1031 colormap = 'cool'
1032 plot_type = 'pcolorbuffer'
1033
1034 def setup(self):
1035 self.xaxis = 'time'
1036 self.ncols = 1
1037 self.nrows = self.data.shape('outlier')[0]
1038 self.nplots = self.nrows
1039 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1040
1041
1042 if not self.xlabel:
1043 self.xlabel = 'Time'
1044
1045 self.ylabel = 'Height [km]'
1046 if not self.titles:
1047 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1048
1049 def update(self, dataOut):
1050
1051 data = {}
1052 data['outlier'] = dataOut.data_outlier
1053
1054 meta = {}
1055
1056 return data, meta
1057
1058 def plot(self):
1059 # self.data.normalize_heights()
1060 self.x = self.data.times
1061 self.y = self.data.yrange
1062 self.z = self.data['outlier']
1063
1064 #self.z = numpy.ma.masked_invalid(self.z)
1065
1066 if self.decimation is None:
1067 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1068 else:
1069 x, y, z = self.fill_gaps(*self.decimate())
1070
1071 for n, ax in enumerate(self.axes):
1072
1073 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1074 self.z[n])
1075 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1076 self.z[n])
1077 data = self.data[-1]
1078 if ax.firsttime:
1079 if self.zlimits is not None:
1080 self.zmin, self.zmax = self.zlimits[n]
1081
1082 ax.plt = ax.pcolormesh(x, y, z[n].T,
1083 vmin=self.zmin,
1084 vmax=self.zmax,
1085 cmap=self.cmaps[n]
1086 )
1087 if self.showprofile:
1088 ax.plot_profile = self.pf_axes[n].plot(data['outlier'][n], self.y)[0]
1089 self.pf_axes[n].set_xlabel('')
1090 else:
1091 if self.zlimits is not None:
1092 self.zmin, self.zmax = self.zlimits[n]
1093 ax.collections.remove(ax.collections[0])
1094 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1095 vmin=self.zmin,
1096 vmax=self.zmax,
1097 cmap=self.cmaps[n]
1098 )
1099 if self.showprofile:
1100 ax.plot_profile.set_data(data['outlier'][n], self.y)
1101 self.pf_axes[n].set_xlabel('')
1102
1103 # class NoiseRTIPlot(Plot):
1104 # '''
1105 # Plot for data_xxxx object
1106 # '''
1107 #
1108 # CODE = 'noise'
1109 # colormap = 'jet'
1110 # plot_type = 'pcolorbuffer'
1111 #
1112 # def setup(self):
1113 # self.xaxis = 'time'
1114 # self.ncols = 1
1115 # self.nrows = self.data.shape('noise')[0]
1116 # self.nplots = self.nrows
1117 # self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1118 #
1119 #
1120 # if not self.xlabel:
1121 # self.xlabel = 'Time'
1122 #
1123 # self.ylabel = 'Height [km]'
1124 # if not self.titles:
1125 # self.titles = ['noise Ch:{}'.format(x) for x in range(self.nrows)]
1126 #
1127 # def update(self, dataOut):
1128 #
1129 # data = {}
1130 # norm = dataOut.max_nIncohInt*dataOut.nProfiles* dataOut.nCohInt*dataOut.windowOfFilter
1131 # print("max incoh: ",dataOut.max_nIncohInt )
1132 # n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1133 # data['noise'] = n0
1134 #
1135 # meta = {}
1136 #
1137 # return data, meta
1138 #
1139 # def plot(self):
1140 # # self.data.normalize_heights()
1141 # self.x = self.data.times
1142 # self.y = self.data.yrange
1143 # self.z = self.data['noise']
1144 #
1145 # #self.z = numpy.ma.masked_invalid(self.z)
1146 #
1147 # if self.decimation is None:
1148 # x, y, z = self.fill_gaps(self.x, self.y, self.z)
1149 # else:
1150 # x, y, z = self.fill_gaps(*self.decimate())
1151 #
1152 # for n, ax in enumerate(self.axes):
1153 #
1154 # self.zmax = self.zmax if self.zmax is not None else numpy.max(
1155 # self.z[n])
1156 # self.zmin = self.zmin if self.zmin is not None else numpy.min(
1157 # self.z[n])
1158 # data = self.data[-1]
1159 # if ax.firsttime:
1160 # if self.zlimits is not None:
1161 # self.zmin, self.zmax = self.zlimits[n]
1162 #
1163 # ax.plt = ax.pcolormesh(x, y, z[n].T,
1164 # vmin=self.zmin,
1165 # vmax=self.zmax,
1166 # cmap=self.cmaps[n]
1167 # )
1168 # if self.showprofile:
1169 # ax.plot_profile = self.pf_axes[n].plot(data['noise'][n], self.y)[0]
1170 # else:
1171 # if self.zlimits is not None:
1172 # self.zmin, self.zmax = self.zlimits[n]
1173 # ax.collections.remove(ax.collections[0])
1174 # ax.plt = ax.pcolormesh(x, y, z[n].T ,
1175 # vmin=self.zmin,
1176 # vmax=self.zmax,
1177 # cmap=self.cmaps[n]
1178 # )
1179 # if self.showprofile:
1180 # ax.plot_profile.set_data(data['noise'][n], self.y)
Show file before | Show file after | Hide comments
@@ -185,6 +185,7 class AMISRReader(ProcessingUnit):
185 self.nsa = self.nsamplesPulse[0,0] #ngates
185 self.nsa = self.nsamplesPulse[0,0] #ngates
186 self.nchannels = len(self.beamCode)
186 self.nchannels = len(self.beamCode)
187 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
187 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
188 #print("IPPS secs: ",self.ippSeconds)
188 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
189 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
189 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
190 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
190
191
Show file before | Show file after | Hide comments
@@ -20,7 +20,7 from schainpy.model.data.jrodata import hildebrand_sekhon
20 from schainpy.model.data import _noise
20 from schainpy.model.data import _noise
21
21
22 from schainpy.utils import log
22 from schainpy.utils import log
23
23 import matplotlib.pyplot as plt
24 from scipy.optimize import curve_fit
24 from scipy.optimize import curve_fit
25
25
26 class SpectraProc(ProcessingUnit):
26 class SpectraProc(ProcessingUnit):
@@ -125,7 +125,7 class SpectraProc(ProcessingUnit):
125 self.dataOut.flagShiftFFT = False
125 self.dataOut.flagShiftFFT = False
126
126
127 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
127 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
128
128 #print("run spc proc")
129 if self.dataIn.type == "Spectra":
129 if self.dataIn.type == "Spectra":
130
130
131 try:
131 try:
@@ -199,6 +199,7 class SpectraProc(ProcessingUnit):
199 if self.profIndex == nProfiles:
199 if self.profIndex == nProfiles:
200
200
201 self.__updateSpecFromVoltage()
201 self.__updateSpecFromVoltage()
202
202 if pairsList == None:
203 if pairsList == None:
203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 else:
205 else:
@@ -208,6 +209,7 class SpectraProc(ProcessingUnit):
208 self.firstdatatime = None
209 self.firstdatatime = None
209 self.profIndex = 0
210 self.profIndex = 0
210
211
212
211 else:
213 else:
212 raise ValueError("The type of input object '%s' is not valid".format(
214 raise ValueError("The type of input object '%s' is not valid".format(
213 self.dataIn.type))
215 self.dataIn.type))
@@ -529,55 +531,56 class getNoise(Operation):
529 # validacion de velocidades
531 # validacion de velocidades
530 indminPoint = None
532 indminPoint = None
531 indmaxPoint = None
533 indmaxPoint = None
534 if self.dataOut.type == 'Spectra':
535 if minVel == None and maxVel == None :
532
536
533 if minVel == None and maxVel == None:
537 freqrange = self.dataOut.getFreqRange(1)
534
538
535 freqrange = self.dataOut.getFreqRange(1)
539 if minFreq == None:
540 minFreq = freqrange[0]
536
541
537 if minFreq == None:
542 if maxFreq == None:
538 minFreq = freqrange[0]
543 maxFreq = freqrange[-1]
539
544
540 if maxFreq == None:
545 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
541 maxFreq = freqrange[-1]
546 if self.warnings:
547 print('minFreq: %.2f is out of the frequency range' % (minFreq))
548 print('minFreq is setting to %.2f' % (freqrange[0]))
549 minFreq = freqrange[0]
542
550
543 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
551 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
544 if self.warnings:
552 if self.warnings:
545 print('minFreq: %.2f is out of the frequency range' % (minFreq))
553 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
546 print('minFreq is setting to %.2f' % (freqrange[0]))
554 print('maxFreq is setting to %.2f' % (freqrange[-1]))
547 minFreq = freqrange[0]
555 maxFreq = freqrange[-1]
548
556
549 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
557 indminPoint = numpy.where(freqrange >= minFreq)
550 if self.warnings:
558 indmaxPoint = numpy.where(freqrange <= maxFreq)
551 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
552 print('maxFreq is setting to %.2f' % (freqrange[-1]))
553 maxFreq = freqrange[-1]
554
559
555 indminPoint = numpy.where(freqrange >= minFreq)
560 else:
556 indmaxPoint = numpy.where(freqrange <= maxFreq)
557
561
558 else:
562 velrange = self.dataOut.getVelRange(1)
559 velrange = self.dataOut.getVelRange(1)
560
563
561 if minVel == None:
564 if minVel == None:
562 minVel = velrange[0]
565 minVel = velrange[0]
563
566
564 if maxVel == None:
567 if maxVel == None:
565 maxVel = velrange[-1]
568 maxVel = velrange[-1]
566
569
567 if (minVel < velrange[0]) or (minVel > maxVel):
570 if (minVel < velrange[0]) or (minVel > maxVel):
568 if self.warnings:
571 if self.warnings:
569 print('minVel: %.2f is out of the velocity range' % (minVel))
572 print('minVel: %.2f is out of the velocity range' % (minVel))
570 print('minVel is setting to %.2f' % (velrange[0]))
573 print('minVel is setting to %.2f' % (velrange[0]))
571 minVel = velrange[0]
574 minVel = velrange[0]
572
575
573 if (maxVel > velrange[-1]) or (maxVel < minVel):
576 if (maxVel > velrange[-1]) or (maxVel < minVel):
574 if self.warnings:
577 if self.warnings:
575 print('maxVel: %.2f is out of the velocity range' % (maxVel))
578 print('maxVel: %.2f is out of the velocity range' % (maxVel))
576 print('maxVel is setting to %.2f' % (velrange[-1]))
579 print('maxVel is setting to %.2f' % (velrange[-1]))
577 maxVel = velrange[-1]
580 maxVel = velrange[-1]
578
581
579 indminPoint = numpy.where(velrange >= minVel)
582 indminPoint = numpy.where(velrange >= minVel)
580 indmaxPoint = numpy.where(velrange <= maxVel)
583 indmaxPoint = numpy.where(velrange <= maxVel)
581
584
582
585
583 # seleccion de indices para rango
586 # seleccion de indices para rango
@@ -606,23 +609,30 class getNoise(Operation):
606 maxIndex = self.dataOut.nHeights - 1
609 maxIndex = self.dataOut.nHeights - 1
607 #############################################################3
610 #############################################################3
608 # seleccion de indices para velocidades
611 # seleccion de indices para velocidades
612 if self.dataOut.type == 'Spectra':
613 try:
614 minIndexFFT = indminPoint[0][0]
615 except:
616 minIndexFFT = 0
609
617
610 try:
618 try:
611 minIndexFFT = indminPoint[0][0]
619 maxIndexFFT = indmaxPoint[0][-1]
612 except:
620 except:
613 minIndexFFT = 0
621 maxIndexFFT = len( self.dataOut.getFreqRange(1))
614
622
615 try:
616 maxIndexFFT = indmaxPoint[0][-1]
617 except:
618 maxIndexFFT = len( self.dataOut.getFreqRange(1))
619
623
620 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
621 self.dataOut.noise_estimation = None
624 self.dataOut.noise_estimation = None
622 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
625 noise = None
623
626 if self.dataOut.type == 'Voltage':
627 noise = self.dataOut.getNoise(ymin_index=minIndex, ymax_index=maxIndex)
628 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
629 elif self.dataOut.type == 'Spectra':
630 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
631 else:
632 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
624 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
633 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
625 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
634 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
635
626 return self.dataOut
636 return self.dataOut
627
637
628
638
@@ -705,7 +715,7 class CleanRayleigh(Operation):
705
715
706
716
707 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
717 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
708
718 #print("runing cleanRayleigh")
709 if not self.isConfig :
719 if not self.isConfig :
710
720
711 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
721 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
@@ -736,7 +746,10 class CleanRayleigh(Operation):
736 if numpy.any(jspc) :
746 if numpy.any(jspc) :
737 #print( jspc.shape, jcspc.shape)
747 #print( jspc.shape, jcspc.shape)
738 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
748 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
739 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
749 try:
750 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
751 except:
752 print("no cspc")
740 self.__dataReady = False
753 self.__dataReady = False
741 #print( jspc.shape, jcspc.shape)
754 #print( jspc.shape, jcspc.shape)
742 dataOut.flagNoData = False
755 dataOut.flagNoData = False
@@ -748,8 +761,11 class CleanRayleigh(Operation):
748 #print( len(self.buffer))
761 #print( len(self.buffer))
749 if numpy.any(self.buffer):
762 if numpy.any(self.buffer):
750 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
763 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
751 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
764 try:
752 self.buffer3 += dataOut.data_dc
765 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
766 self.buffer3 += dataOut.data_dc
767 except:
768 pass
753 else:
769 else:
754 self.buffer = dataOut.data_spc
770 self.buffer = dataOut.data_spc
755 self.buffer2 = dataOut.data_cspc
771 self.buffer2 = dataOut.data_cspc
@@ -762,7 +778,9 class CleanRayleigh(Operation):
762
778
763
779
764 #index = tini.tm_hour*12+tini.tm_min/5
780 #index = tini.tm_hour*12+tini.tm_min/5
765 '''REVISAR'''
781 '''
782 REVISAR
783 '''
766 # jspc = jspc/self.nFFTPoints/self.normFactor
784 # jspc = jspc/self.nFFTPoints/self.normFactor
767 # jcspc = jcspc/self.nFFTPoints/self.normFactor
785 # jcspc = jcspc/self.nFFTPoints/self.normFactor
768
786
@@ -776,6 +794,7 class CleanRayleigh(Operation):
776
794
777 dataOut.data_dc = self.buffer3
795 dataOut.data_dc = self.buffer3
778 dataOut.nIncohInt *= self.nIntProfiles
796 dataOut.nIncohInt *= self.nIntProfiles
797 dataOut.max_nIncohInt = self.nIntProfiles
779 dataOut.utctime = self.currentTime #tiempo promediado
798 dataOut.utctime = self.currentTime #tiempo promediado
780 #print("Time: ",time.localtime(dataOut.utctime))
799 #print("Time: ",time.localtime(dataOut.utctime))
781 # dataOut.data_spc = sat_spectra
800 # dataOut.data_spc = sat_spectra
@@ -787,7 +806,7 class CleanRayleigh(Operation):
787 return dataOut
806 return dataOut
788
807
789 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
808 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
790 #print("OP cleanRayleigh")
809 print("OP cleanRayleigh")
791 #import matplotlib.pyplot as plt
810 #import matplotlib.pyplot as plt
792 #for k in range(149):
811 #for k in range(149):
793 #channelsProcssd = []
812 #channelsProcssd = []
@@ -803,6 +822,8 class CleanRayleigh(Operation):
803
822
804 ###ONLY FOR TEST:
823 ###ONLY FOR TEST:
805 raxs = math.ceil(math.sqrt(self.nPairs))
824 raxs = math.ceil(math.sqrt(self.nPairs))
825 if raxs == 0:
826 raxs = 1
806 caxs = math.ceil(self.nPairs/raxs)
827 caxs = math.ceil(self.nPairs/raxs)
807 if self.nPairs <4:
828 if self.nPairs <4:
808 raxs = 2
829 raxs = 2
@@ -1100,12 +1121,13 class IntegrationFaradaySpectra(Operation):
1100 __dataReady = False
1121 __dataReady = False
1101
1122
1102 __timeInterval = None
1123 __timeInterval = None
1103
1124 n_ints = None #matriz de numero de integracions (CH,HEI)
1104 n = None
1125 n = None
1105 minHei_ind = None
1126 minHei_ind = None
1106 maxHei_ind = None
1127 maxHei_ind = None
1107 navg = 1.0
1128 navg = 1.0
1108 factor = 0.0
1129 factor = 0.0
1130 dataoutliers = None # (CHANNELS, HEIGHTS)
1109
1131
1110 def __init__(self):
1132 def __init__(self):
1111
1133
@@ -1138,6 +1160,7 class IntegrationFaradaySpectra(Operation):
1138 self.navg = avg
1160 self.navg = avg
1139 #self.ByLags = dataOut.ByLags ###REDEFINIR
1161 #self.ByLags = dataOut.ByLags ###REDEFINIR
1140 self.ByLags = False
1162 self.ByLags = False
1163 self.maxProfilesInt = 1
1141
1164
1142 if DPL != None:
1165 if DPL != None:
1143 self.DPL=DPL
1166 self.DPL=DPL
@@ -1251,6 +1274,8 class IntegrationFaradaySpectra(Operation):
1251 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1274 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1252 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1275 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1253
1276
1277 self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
1278
1254 for k in range(self.minHei_ind,self.maxHei_ind):
1279 for k in range(self.minHei_ind,self.maxHei_ind):
1255 try:
1280 try:
1256 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1281 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
@@ -1266,22 +1291,31 class IntegrationFaradaySpectra(Operation):
1266 #sortIDs=[]
1291 #sortIDs=[]
1267 outliers_IDs=[]
1292 outliers_IDs=[]
1268
1293
1269 for j in range(self.nProfiles):
1294 for j in range(self.nProfiles): #frecuencias en el tiempo
1270 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1295 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1271 # continue
1296 # continue
1272 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1297 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1273 # continue
1298 # continue
1274 buffer=buffer1[:,j]
1299 buffer=buffer1[:,j]
1275 sortdata = numpy.sort(buffer, axis=None)
1300 sortdata = numpy.sort(buffer, axis=None)
1301
1276 sortID=buffer.argsort()
1302 sortID=buffer.argsort()
1277 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1303 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1278
1304
1279 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1305 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1280
1306
1307 # fig,ax = plt.subplots()
1308 # ax.set_title(str(k)+" "+str(j))
1309 # x=range(len(sortdata))
1310 # ax.scatter(x,sortdata)
1311 # ax.axvline(index)
1312 # plt.show()
1313
1281 indexes.append(index)
1314 indexes.append(index)
1282 #sortIDs.append(sortID)
1315 #sortIDs.append(sortID)
1283 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1316 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1284
1317
1318 #print("Outliers: ",outliers_IDs)
1285 outliers_IDs=numpy.array(outliers_IDs)
1319 outliers_IDs=numpy.array(outliers_IDs)
1286 outliers_IDs=outliers_IDs.ravel()
1320 outliers_IDs=outliers_IDs.ravel()
1287 outliers_IDs=numpy.unique(outliers_IDs)
1321 outliers_IDs=numpy.unique(outliers_IDs)
@@ -1289,31 +1323,45 class IntegrationFaradaySpectra(Operation):
1289 indexes=numpy.array(indexes)
1323 indexes=numpy.array(indexes)
1290 indexmin=numpy.min(indexes)
1324 indexmin=numpy.min(indexes)
1291
1325
1326
1327 #print(indexmin,buffer1.shape[0], k)
1328
1329 # fig,ax = plt.subplots()
1330 # ax.plot(sortdata)
1331 # ax2 = ax.twinx()
1332 # x=range(len(indexes))
1333 # #plt.scatter(x,indexes)
1334 # ax2.scatter(x,indexes)
1335 # plt.show()
1336
1292 if indexmin != buffer1.shape[0]:
1337 if indexmin != buffer1.shape[0]:
1293 if self.nChannels > 1:
1338 if self.nChannels > 1:
1294 cspc_outliers_exist= True
1339 cspc_outliers_exist= True
1295
1340
1296 lt=outliers_IDs
1341 lt=outliers_IDs
1297 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1342 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1298
1343
1299 for p in list(outliers_IDs):
1344 for p in list(outliers_IDs):
1300 buffer1[p,:]=avg
1345 #buffer1[p,:]=avg
1346 buffer1[p,:] = numpy.NaN
1347
1348 self.dataOutliers[i,k] = len(outliers_IDs)
1301
1349
1302 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1350 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1303 ###cspc IDs
1351
1304 #indexmin_cspc+=indexmin_cspc
1305 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1352 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1306
1353
1307 #if not breakFlag:
1354
1355
1308 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1356 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1309 if cspc_outliers_exist :
1357 if cspc_outliers_exist :
1310 #sortdata=numpy.sort(buffer_cspc,axis=0)
1358
1311 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1312 lt=outliers_IDs_cspc
1359 lt=outliers_IDs_cspc
1313
1360
1314 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1361 #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1315 for p in list(outliers_IDs_cspc):
1362 for p in list(outliers_IDs_cspc):
1316 buffer_cspc[p,:]=avg
1363 #buffer_cspc[p,:]=avg
1364 buffer_cspc[p,:] = numpy.NaN
1317
1365
1318 try:
1366 try:
1319 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1367 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
@@ -1325,39 +1373,36 class IntegrationFaradaySpectra(Operation):
1325
1373
1326
1374
1327
1375
1328
1376 nOutliers = len(outliers_IDs)
1377 #print("Outliers n: ",self.dataOutliers,nOutliers)
1329 buffer=None
1378 buffer=None
1330 bufferH=None
1379 bufferH=None
1331 buffer1=None
1380 buffer1=None
1332 buffer_cspc=None
1381 buffer_cspc=None
1333
1382
1334 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1335 #exit()
1336
1383
1337 buffer=None
1384 buffer=None
1338 #print(self.__buffer_spc[:,1,3,20,0])
1385
1339 #print(self.__buffer_spc[:,1,5,37,0])
1386 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1340 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1387 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1341 try:
1388 try:
1342 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1389 #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1390 data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
1343 except:
1391 except:
1344 #print("No cpsc",e)
1392 #print("No cpsc",e)
1345 pass
1393 pass
1346
1394
1347
1395
1348 #print(numpy.shape(data_spc))
1349 #data_spc[1,4,20,0]=numpy.nan
1350
1351 #data_cspc = self.__buffer_cspc
1352 #print("pushData pre Done")
1353 data_dc = self.__buffer_dc
1396 data_dc = self.__buffer_dc
1354 n = self.__profIndex
1397 #(CH, HEIGH)
1398 self.maxProfilesInt = self.__profIndex
1399 n = self.__profIndex - self.dataOutliers
1355
1400
1356 self.__buffer_spc = []
1401 self.__buffer_spc = []
1357 self.__buffer_cspc = []
1402 self.__buffer_cspc = []
1358 self.__buffer_dc = 0
1403 self.__buffer_dc = 0
1359 self.__profIndex = 0
1404 self.__profIndex = 0
1360 #print("pushData Done")
1405
1361 return data_spc, data_cspc, data_dc, n
1406 return data_spc, data_cspc, data_dc, n
1362
1407
1363 def byProfiles(self, *args):
1408 def byProfiles(self, *args):
@@ -1372,7 +1417,7 class IntegrationFaradaySpectra(Operation):
1372 if self.__profIndex == self.n:
1417 if self.__profIndex == self.n:
1373
1418
1374 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1419 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1375 self.n = n
1420 self.n_ints = n
1376 self.__dataReady = True
1421 self.__dataReady = True
1377
1422
1378 return avgdata_spc, avgdata_cspc, avgdata_dc
1423 return avgdata_spc, avgdata_cspc, avgdata_dc
@@ -1388,7 +1433,7 class IntegrationFaradaySpectra(Operation):
1388
1433
1389 if (datatime - self.__initime) >= self.__integrationtime:
1434 if (datatime - self.__initime) >= self.__integrationtime:
1390 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1435 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1391 self.n = n
1436 self.n_ints = n
1392 self.__dataReady = True
1437 self.__dataReady = True
1393
1438
1394 return avgdata_spc, avgdata_cspc, avgdata_dc
1439 return avgdata_spc, avgdata_cspc, avgdata_dc
@@ -1450,7 +1495,7 class IntegrationFaradaySpectra(Operation):
1450 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1495 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1451 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1496 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1452 self.dataOut.data_dc = avgdata_dc
1497 self.dataOut.data_dc = avgdata_dc
1453
1498 self.dataOut.data_outlier = self.dataOutliers
1454
1499
1455 else:
1500 else:
1456 self.dataOut.dataLag_spc = avgdata_spc
1501 self.dataOut.dataLag_spc = avgdata_spc
@@ -1462,10 +1507,12 class IntegrationFaradaySpectra(Operation):
1462 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1507 self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
1463
1508
1464
1509
1465 self.dataOut.nIncohInt *= self.n
1510 self.dataOut.nIncohInt *= self.n_ints
1511 self.dataOut.max_nIncohInt = self.maxProfilesInt
1512 #print(self.dataOut.max_nIncohInt)
1466 self.dataOut.utctime = avgdatatime
1513 self.dataOut.utctime = avgdatatime
1467 self.dataOut.flagNoData = False
1514 self.dataOut.flagNoData = False
1468
1515 #print("Faraday Integration DONE...")
1469 return self.dataOut
1516 return self.dataOut
1470
1517
1471 class removeInterference(Operation):
1518 class removeInterference(Operation):
@@ -1705,7 +1752,8 class IncohInt(Operation):
1705 __dataReady = False
1752 __dataReady = False
1706
1753
1707 __timeInterval = None
1754 __timeInterval = None
1708
1755 incohInt = 0
1756 nOutliers = 0
1709 n = None
1757 n = None
1710
1758
1711 def __init__(self):
1759 def __init__(self):
@@ -1734,7 +1782,8 class IncohInt(Operation):
1734 self.__profIndex = 0
1782 self.__profIndex = 0
1735 self.__dataReady = False
1783 self.__dataReady = False
1736 self.__byTime = False
1784 self.__byTime = False
1737
1785 self.incohInt = 0
1786 self.nOutliers = 0
1738 if n is None and timeInterval is None:
1787 if n is None and timeInterval is None:
1739 raise ValueError("n or timeInterval should be specified ...")
1788 raise ValueError("n or timeInterval should be specified ...")
1740
1789
@@ -1788,7 +1837,7 class IncohInt(Operation):
1788 self.__buffer_spc = 0
1837 self.__buffer_spc = 0
1789 self.__buffer_cspc = 0
1838 self.__buffer_cspc = 0
1790 self.__buffer_dc = 0
1839 self.__buffer_dc = 0
1791 self.__profIndex = 0
1840
1792
1841
1793 return data_spc, data_cspc, data_dc, n
1842 return data_spc, data_cspc, data_dc, n
1794
1843
@@ -1845,6 +1894,9 class IncohInt(Operation):
1845 if n == 1:
1894 if n == 1:
1846 return dataOut
1895 return dataOut
1847
1896
1897 if dataOut.flagNoData == True:
1898 return dataOut
1899
1848 dataOut.flagNoData = True
1900 dataOut.flagNoData = True
1849
1901
1850 if not self.isConfig:
1902 if not self.isConfig:
@@ -1855,15 +1907,21 class IncohInt(Operation):
1855 dataOut.data_spc,
1907 dataOut.data_spc,
1856 dataOut.data_cspc,
1908 dataOut.data_cspc,
1857 dataOut.data_dc)
1909 dataOut.data_dc)
1858
1910 self.incohInt += dataOut.nIncohInt
1911 self.nOutliers += dataOut.data_outlier
1859 if self.__dataReady:
1912 if self.__dataReady:
1860
1913
1861 dataOut.data_spc = avgdata_spc
1914 dataOut.data_spc = avgdata_spc
1862 dataOut.data_cspc = avgdata_cspc
1915 dataOut.data_cspc = avgdata_cspc
1863 dataOut.data_dc = avgdata_dc
1916 dataOut.data_dc = avgdata_dc
1864 dataOut.nIncohInt *= self.n
1917 dataOut.nIncohInt = self.incohInt
1918 dataOut.data_outlier = self.nOutliers
1865 dataOut.utctime = avgdatatime
1919 dataOut.utctime = avgdatatime
1866 dataOut.flagNoData = False
1920 dataOut.flagNoData = False
1921 dataOut.max_nIncohInt = self.__profIndex
1922 self.incohInt = 0
1923 self.nOutliers = 0
1924 self.__profIndex = 0
1867
1925
1868 return dataOut
1926 return dataOut
1869
1927
Show file before | Show file after | Hide comments
@@ -8,7 +8,8 from schainpy.model.io.utils import getHei_index
8 from time import time
8 from time import time
9 import datetime
9 import datetime
10 import numpy
10 import numpy
11 import copy
11 #import copy
12 from schainpy.model.data import _noise
12
13
13 class VoltageProc(ProcessingUnit):
14 class VoltageProc(ProcessingUnit):
14
15
@@ -1759,6 +1760,8 class SSheightProfiles2(Operation):
1759 profileShape = None
1760 profileShape = None
1760 sshProfiles = None
1761 sshProfiles = None
1761 profileIndex = None
1762 profileIndex = None
1763 deltaHeight = None
1764 init_range = None
1762
1765
1763 def __init__(self, **kwargs):
1766 def __init__(self, **kwargs):
1764
1767
@@ -1780,7 +1783,8 class SSheightProfiles2(Operation):
1780 if residue != 0:
1783 if residue != 0:
1781 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))
1784 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))
1782
1785
1783 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1786 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1787 self.init_range = dataOut.heightList[0]
1784 #numberProfile = self.nsamples
1788 #numberProfile = self.nsamples
1785 self.numberSamples = (self.__nHeis - self.nsamples)/self.step
1789 self.numberSamples = (self.__nHeis - self.nsamples)/self.step
1786
1790
@@ -1800,21 +1804,25 class SSheightProfiles2(Operation):
1800
1804
1801 if repeat is not None:
1805 if repeat is not None:
1802 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1806 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1803 #print("buff, data, :",self.buffer.shape, data.shape)
1807 if data.ndim == 2:
1808 data = data.reshape(1,1,self.__nHeis )
1809 #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape)
1804 for i in range(int(self.new_nHeights)): #nuevas alturas
1810 for i in range(int(self.new_nHeights)): #nuevas alturas
1805 if code is not None:
1811 if code is not None:
1806 self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]*code_block
1812 self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]*code_block
1807 else:
1813 else:
1808 self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1814 self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1809
1815
1810 for j in range(self.__nChannels): #en los cananles
1816 for j in range(self.__nChannels): #en los cananles
1811 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1817 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
1812
1818 #print("new profs Done")
1813
1819
1814
1820
1815
1821
1816 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1822 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1817
1823
1824 if dataOut.flagNoData == True:
1825 return dataOut
1818 dataOut.flagNoData = True
1826 dataOut.flagNoData = True
1819 #print("init data shape:", dataOut.data.shape)
1827 #print("init data shape:", dataOut.data.shape)
1820 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
1828 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
@@ -1838,24 +1846,26 class SSheightProfiles2(Operation):
1838 #print("dataOut nProfiles:", dataOut.nProfiles)
1846 #print("dataOut nProfiles:", dataOut.nProfiles)
1839 for profile in range(nprof):
1847 for profile in range(nprof):
1840 if dataOut.flagDataAsBlock:
1848 if dataOut.flagDataAsBlock:
1849 #print("read blocks")
1841 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
1850 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
1842 else:
1851 else:
1843 self.getNewProfiles(dataOut.data, code=code, repeat=repeat)
1852 #print("read profiles")
1853 self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe
1844 if profile == 0:
1854 if profile == 0:
1845 dataBlock = self.sshProfiles.copy()
1855 dataBlock = self.sshProfiles.copy()
1846 else: #by blocks
1856 else: #by blocks
1847 dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
1857 dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
1848 print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
1858 #print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
1849
1859
1850 profileIndex = self.nsamples
1860 profileIndex = self.nsamples
1851 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1861 #deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1852 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1862 ippSeconds = (self.deltaHeight*1.0e-6)/(0.15)
1853
1863
1854
1864
1855 dataOut.data = dataBlock
1865 dataOut.data = dataBlock
1856 dataOut.heightList = numpy.arange(self.new_nHeights) *self.step*deltaHeight + dataOut.heightList[0]
1866 #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
1857 #print("show me: ",dataOut.nHeights, self.new_nHeights)
1867 dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
1858 #dataOut.nHeights = int(self.new_nHeights)
1868
1859 dataOut.ippSeconds = ippSeconds
1869 dataOut.ippSeconds = ippSeconds
1860 dataOut.step = self.step
1870 dataOut.step = self.step
1861 dataOut.flagNoData = False
1871 dataOut.flagNoData = False
@@ -1877,17 +1887,17 class SSheightProfiles2(Operation):
1877
1887
1878 #import skimage.color
1888 #import skimage.color
1879 #import skimage.io
1889 #import skimage.io
1880 import matplotlib.pyplot as plt
1890 #import matplotlib.pyplot as plt
1881
1891
1882 class clean2DArray(Operation):
1892 class removeProfileByFaradayHS(Operation):
1883 '''
1893 '''
1884
1894
1885 '''
1895 '''
1886 isConfig = False
1896 isConfig = False
1887 n = None
1897 n = None
1888 __timeInterval = None
1898
1889 __profIndex = 0
1899 __profIndex = 0
1890 __byTime = False
1900
1891 __dataReady = False
1901 __dataReady = False
1892 __buffer_data = []
1902 __buffer_data = []
1893 __buffer_times = []
1903 __buffer_times = []
@@ -1901,6 +1911,7 class clean2DArray(Operation):
1901 end_prof = 0
1911 end_prof = 0
1902 n_profiles = 0
1912 n_profiles = 0
1903 first_utcBlock = None
1913 first_utcBlock = None
1914 outliers_IDs_list = []
1904 __dh = 0
1915 __dh = 0
1905
1916
1906 def __init__(self, **kwargs):
1917 def __init__(self, **kwargs):
@@ -1908,42 +1919,122 class clean2DArray(Operation):
1908 Operation.__init__(self, **kwargs)
1919 Operation.__init__(self, **kwargs)
1909 self.isConfig = False
1920 self.isConfig = False
1910
1921
1911 def setup(self,dataOut, n=None , timeInterval=None):
1922 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3, minHei=None, maxHei=None):
1912
1923
1913 if n == None and timeInterval == None:
1924 if n == None and timeInterval == None:
1914 raise ValueError("nprofiles or timeInterval should be specified ...")
1925 raise ValueError("nprofiles or timeInterval should be specified ...")
1915
1926
1916 if n != None:
1927 if n != None:
1917 self.n = n
1928 self.n = n
1918 self.__byTime = False
1919 else:
1920 self.__timeInterval = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
1921 self.n = 9999
1922 self.__byTime = True
1923
1929
1930 self.navg = navg
1931 self.profileMargin = profileMargin
1932 self.thHistOutlier = thHistOutlier
1924 self.__profIndex = 0
1933 self.__profIndex = 0
1925 self.buffer = None
1934 self.buffer = None
1926 self.lenProfileOut = 1
1935 self.lenProfileOut = 1
1927
1936 self.n_prof_released = 0
1937 self.heightList = dataOut.heightList
1928 self.init_prof = 0
1938 self.init_prof = 0
1929 self.end_prof = 0
1939 self.end_prof = 0
1930 self.n_profiles = 0
1940 self.n_profiles = 0
1931 self.first_utcBlock = None
1941 self.first_utcBlock = None
1932 self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
1942 self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
1943 minHei = minHei
1944 maxHei = maxHei
1945 if minHei==None :
1946 minHei = dataOut.heightList[0]
1947 if maxHei==None :
1948 maxHei = dataOut.heightList[-1]
1949 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
1950
1951
1952 def filterSatsProfiles(self):
1953 data = self.__buffer_data
1954 #print(data.shape)
1955 nChannels, profiles, heights = data.shape
1956 indexes=[]
1957 outliers_IDs=[]
1958 for c in range(nChannels):
1959 for h in range(self.minHei_idx, self.maxHei_idx):
1960 power = data[c,:,h] * numpy.conjugate(data[c,:,h])
1961 power = power.real
1962 #power = (numpy.abs(data[c,:,h].real))
1963 sortdata = numpy.sort(power, axis=None)
1964 sortID=power.argsort()
1965 index = _noise.hildebrand_sekhon2(sortdata,self.navg) #0.75-> buen valor
1966
1967 indexes.append(index)
1968 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1969 # print(outliers_IDs)
1970 # fig,ax = plt.subplots()
1971 # #ax.set_title(str(k)+" "+str(j))
1972 # x=range(len(sortdata))
1973 # ax.scatter(x,sortdata)
1974 # ax.axvline(index)
1975 # plt.grid()
1976 # plt.show()
1977
1978 outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
1979 outliers_IDs = numpy.unique(outliers_IDs)
1980 outs_lines = numpy.sort(outliers_IDs)
1981 # #print("outliers Ids: ", outs_lines, outs_lines.shape)
1982 #hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
1983
1984
1985 #Agrupando el histograma de outliers,
1986 my_bins = numpy.linspace(0,9600, 96, endpoint=False)
1987
1988
1989 hist, bins = numpy.histogram(outs_lines,bins=my_bins)
1990 hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier
1991 #print(hist_outliers_indexes[0])
1992 bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
1993 #print(bins_outliers_indexes)
1994 outlier_loc_index = []
1995
1996 #outlier_loc_index = [k for k in range(bins_outliers_indexes[n]-50,bins_outliers_indexes[n+1]+50) for n in range(len(bins_outliers_indexes)-1) ]
1997 for n in range(len(bins_outliers_indexes)-1):
1998 #outlier_loc_index = [k for k in range(bins_outliers_indexes[n]-50,bins_outliers_indexes[n+1]+50)]
1999 for k in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin):
2000 outlier_loc_index.append(k)
2001
2002 outlier_loc_index = numpy.asarray(outlier_loc_index)
2003 #print(numpy.unique(outlier_loc_index))
2004
2005
2006
2007 # x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
2008 # fig, ax = plt.subplots(1,2,figsize=(8, 6))
2009 #
2010 # dat = data[0,:,:].real
2011 # m = numpy.nanmean(dat)
2012 # o = numpy.nanstd(dat)
2013 # #print(m, o, x.shape, y.shape)
2014 # c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2015 # ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
2016 # fig.colorbar(c)
2017 # ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
2018 # ax[1].hist(outs_lines,bins=my_bins)
2019 # plt.show()
2020
2021
2022 self.outliers_IDs_list = numpy.unique(outlier_loc_index)
2023 return data
1933
2024
1934 def cleanOutliersByBlock(self):
2025 def cleanOutliersByBlock(self):
1935 #print(self.__buffer_data[0].shape)
2026 #print(self.__buffer_data[0].shape)
1936 data = self.__buffer_data#.copy()
2027 data = self.__buffer_data#.copy()
1937 #print("cleaning shape inpt: ",data.shape)
2028 #print("cleaning shape inpt: ",data.shape)
1938
2029 '''
1939 self.__buffer_data = []
2030 self.__buffer_data = []
1940
2031
1941 spectrum = numpy.fft.fft2(data, axes=(0,2))
2032 spectrum = numpy.fft.fft2(data, axes=(0,2))
1942 #print("spc : ",spectrum.shape)
2033 #print("spc : ",spectrum.shape)
1943 (nch,nsamples, nh) = spectrum.shape
2034 (nch,nsamples, nh) = spectrum.shape
1944
2035 data2 = None
1945 #nch =
2036 #print(data.shape)
1946
2037 spectrum2 = spectrum.copy()
1947 for ch in range(nch):
2038 for ch in range(nch):
1948 dh = self.__dh
2039 dh = self.__dh
1949 dt1 = (dh*1.0e-6)/(0.15)
2040 dt1 = (dh*1.0e-6)/(0.15)
@@ -1951,57 +2042,48 class clean2DArray(Operation):
1951 freqv = numpy.fft.fftfreq(nh, d=dt1)
2042 freqv = numpy.fft.fftfreq(nh, d=dt1)
1952 freqh = numpy.fft.fftfreq(self.n, d=dt2)
2043 freqh = numpy.fft.fftfreq(self.n, d=dt2)
1953 #print("spc loop: ")
2044 #print("spc loop: ")
1954 x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
1955 z = numpy.abs(spectrum[ch,:,:])
1956
1957
1958 dat = numpy.log10(z.T)
1959 m = numpy.mean(dat)
1960 o = numpy.std(dat)
1961 fig, ax = plt.subplots(figsize=(8, 6))
1962
2045
1963 c = ax.pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
1964 #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
1965 date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
1966 #strftime('%Y-%m-%d %H:%M:%S')
1967 ax.set_title('Spectrum magnitude '+date_time)
1968 fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
1969 fig.colorbar(c)
1970
2046
1971 plt.show()
1972
2047
2048 x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
2049 z = numpy.abs(spectrum[ch,:,:])
2050 # Find all peaks higher than the 98th percentile
2051 peaks = z < numpy.percentile(z, 98)
2052 #print(peaks)
2053 # Set those peak coefficients to zero
2054 spectrum2 = spectrum2 * peaks.astype(int)
2055 data2 = numpy.fft.ifft2(spectrum2)
1973
2056
2057 dat = numpy.log10(z.T)
2058 dat2 = numpy.log10(spectrum2.T)
2059
2060 # m = numpy.mean(dat)
2061 # o = numpy.std(dat)
2062 # fig, ax = plt.subplots(2,1,figsize=(8, 6))
2063 #
2064 # c = ax[0].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2065 # #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
2066 # date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
2067 # #strftime('%Y-%m-%d %H:%M:%S')
2068 # ax[0].set_title('Spectrum magnitude '+date_time)
2069 # fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
2070 #
2071 #
2072 # c = ax[1].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
2073 # fig.colorbar(c)
2074 # plt.show()
2075
2076 #print(data2.shape)
2077
2078 data = data2
1974
2079
1975 #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
2080 #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
2081 '''
2082 #print("cleanOutliersByBlock Done")
1976
2083
1977 print("cleanOutliersByBlock Done")
2084 return self.filterSatsProfiles()
1978 return data
1979
1980 def byTime(self, data, datatime):
1981
1982 self.__dataReady = False
1983
1984 self.fillBuffer(data, datatime)
1985 dataBlock = None
1986
1987 if (datatime - self.__initime) >= self.__timeInterval:
1988 dataBlock = self.cleanOutliersByBlock()
1989 self.__dataReady = True
1990 self.n = self.__profIndex
1991 return dataBlock
1992
1993 def byProfiles(self, data, datatime):
1994 self.__dataReady = False
1995
1996 self.fillBuffer(data, datatime)
1997 dataBlock = None
1998
2085
1999 if self.__profIndex == self.n:
2000 #print("apnd : ",data)
2001 dataBlock = self.cleanOutliersByBlock()
2002 self.__dataReady = True
2003
2086
2004 return dataBlock
2005
2087
2006 def fillBuffer(self, data, datatime):
2088 def fillBuffer(self, data, datatime):
2007
2089
@@ -2011,54 +2093,60 class clean2DArray(Operation):
2011 else:
2093 else:
2012 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2094 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2013 self.__profIndex += 1
2095 self.__profIndex += 1
2014 self.__buffer_times.append(datatime)
2096 #self.__buffer_times.append(datatime)
2015
2097
2016 def getData(self, data, datatime=None):
2098 def getData(self, data, datatime=None):
2017
2099
2018 if self.__profIndex == 0:
2100 if self.__profIndex == 0:
2019 self.__initime = datatime
2101 self.__initime = datatime
2020
2102
2021 if self.__byTime:
2022 dataBlock = self.byTime(data, datatime)
2023 else:
2024 dataBlock = self.byProfiles(data, datatime)
2025
2103
2104 self.__dataReady = False
2026
2105
2106 self.fillBuffer(data, datatime)
2107 dataBlock = None
2108
2109 if self.__profIndex == self.n:
2110 #print("apnd : ",data)
2111 #dataBlock = self.cleanOutliersByBlock()
2112 dataBlock = self.filterSatsProfiles()
2113 self.__dataReady = True
2114
2115 return dataBlock
2027
2116
2028 if dataBlock is None:
2117 if dataBlock is None:
2029 return None
2118 return None, None
2030
2119
2031
2120
2032
2121
2033 return dataBlock
2122 return dataBlock
2034
2123
2035 def releaseBlock(self, dataOut):
2124 def releaseBlock(self):
2036
2125
2037 if self.n % self.lenProfileOut != 0:
2126 if self.n % self.lenProfileOut != 0:
2038 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2127 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2039 return None
2128 return None
2040
2129
2041 dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2130 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2042 dataOut.profileIndex = self.lenProfileOut
2131
2043 dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds
2044 self.init_prof = self.end_prof
2132 self.init_prof = self.end_prof
2045 self.end_prof += self.lenProfileOut
2133 self.end_prof += self.lenProfileOut
2046 print("data release shape: ",dataOut.data.shape, self.end_prof)
2134 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2135 self.n_prof_released += 1
2047
2136
2048 if self.end_prof >= (self.n +self.lenProfileOut):
2049 self.init_prof = 0
2050 self.__profIndex = 0
2051 self.buffer = None
2052 dataOut.buffer_empty = True
2053 return dataOut
2054
2137
2055 def run(self, dataOut, n=None, timeInterval=None, nProfilesOut=1):
2138 #print("f_no_data ", dataOut.flagNoData)
2139 return data
2140
2141 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None):
2056 #print("run op buffer 2D")
2142 #print("run op buffer 2D")
2057 self.nChannels = dataOut.nChannels
2143 self.nChannels = dataOut.nChannels
2058 self.nHeights = dataOut.nHeights
2144 self.nHeights = dataOut.nHeights
2059
2145
2060 if not self.isConfig:
2146 if not self.isConfig:
2061 self.setup(dataOut,n=n, timeInterval=timeInterval)
2147 #print("init p idx: ", dataOut.profileIndex )
2148 self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,
2149 thHistOutlier=th_hist_outlier,minHei=minHei, maxHei=maxHei)
2062 self.isConfig = True
2150 self.isConfig = True
2063
2151
2064 dataBlock = None
2152 dataBlock = None
@@ -2066,6 +2154,7 class clean2DArray(Operation):
2066 if not dataOut.buffer_empty: #hay datos acumulados
2154 if not dataOut.buffer_empty: #hay datos acumulados
2067
2155
2068 if self.init_prof == 0:
2156 if self.init_prof == 0:
2157 self.n_prof_released = 0
2069 self.lenProfileOut = nProfilesOut
2158 self.lenProfileOut = nProfilesOut
2070 dataOut.flagNoData = False
2159 dataOut.flagNoData = False
2071 #print("tp 2 ",dataOut.data.shape)
2160 #print("tp 2 ",dataOut.data.shape)
@@ -2074,13 +2163,45 class clean2DArray(Operation):
2074 #print("rel: ",self.buffer[:,-1,:])
2163 #print("rel: ",self.buffer[:,-1,:])
2075 self.init_prof = 0
2164 self.init_prof = 0
2076 self.end_prof = self.lenProfileOut
2165 self.end_prof = self.lenProfileOut
2077 self.first_utcBlock = dataOut.utctime
2166
2078 dataOut.nProfiles = self.lenProfileOut
2167 dataOut.nProfiles = self.lenProfileOut
2079 dataOut.error = False
2168 if nProfilesOut == 1:
2169 dataOut.flagDataAsBlock = False
2170 else:
2080 dataOut.flagDataAsBlock = True
2171 dataOut.flagDataAsBlock = True
2081 #print("prof: ",self.init_prof)
2172 #print("prof: ",self.init_prof)
2082 dataOut.flagNoData = False
2173 dataOut.flagNoData = False
2083 return self.releaseBlock(dataOut)
2174 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
2175 #print("omitting: ", self.n_prof_released)
2176 dataOut.flagNoData = True
2177
2178 dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds
2179 #dataOut.data = self.releaseBlock()
2180 #########################################################3
2181 if self.n % self.lenProfileOut != 0:
2182 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2183 return None
2184
2185 dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2186
2187 self.init_prof = self.end_prof
2188 self.end_prof += self.lenProfileOut
2189 #print("data release shape: ",dataOut.data.shape, self.end_prof)
2190 self.n_prof_released += 1
2191
2192 if self.end_prof >= (self.n +self.lenProfileOut):
2193
2194 self.init_prof = 0
2195 self.__profIndex = 0
2196 self.buffer = None
2197 dataOut.buffer_empty = True
2198 self.outliers_IDs_list = []
2199 self.n_prof_released = 0
2200 dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
2201 #print("cleaning...", dataOut.buffer_empty)
2202 dataOut.profileIndex = 0 #self.lenProfileOut
2203 ####################################################################
2204 return dataOut
2084
2205
2085
2206
2086 #print("tp 223 ",dataOut.data.shape)
2207 #print("tp 223 ",dataOut.data.shape)
@@ -2098,13 +2219,16 class clean2DArray(Operation):
2098 sys.exit()
2219 sys.exit()
2099
2220
2100 if self.__dataReady:
2221 if self.__dataReady:
2222 #print("omitting: ", len(self.outliers_IDs_list))
2101 self.__count_exec = 0
2223 self.__count_exec = 0
2102 #dataOut.data =
2224 #dataOut.data =
2103 self.buffer = numpy.flip(dataBlock, axis=1)
2225 #self.buffer = numpy.flip(dataBlock, axis=1)
2104
2226 self.buffer = dataBlock
2227 self.first_utcBlock = self.__initime
2105 dataOut.utctime = self.__initime
2228 dataOut.utctime = self.__initime
2106 dataOut.nProfiles = self.__profIndex
2229 dataOut.nProfiles = self.__profIndex
2107 #dataOut.flagNoData = False
2230 #dataOut.flagNoData = False
2231 self.init_prof = 0
2108 self.__profIndex = 0
2232 self.__profIndex = 0
2109 self.__initime = None
2233 self.__initime = None
2110 dataBlock = None
2234 dataBlock = None
@@ -2112,16 +2236,15 class clean2DArray(Operation):
2112 dataOut.error = False
2236 dataOut.error = False
2113 dataOut.useInputBuffer = True
2237 dataOut.useInputBuffer = True
2114 dataOut.buffer_empty = False
2238 dataOut.buffer_empty = False
2115 print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
2239 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
2116 int(dataOut.nProfiles),int(dataOut.nHeights)))
2240
2117 #return None
2118
2241
2119
2242
2120 #print(self.__count_exec)
2243 #print(self.__count_exec)
2121
2244
2122 return dataOut
2245 return dataOut
2123
2246
2124 class CleanProfileSats(Operation):
2247 class RemoveProfileSats(Operation):
2125 '''
2248 '''
2126 Omite los perfiles contaminados con señal de satelites,
2249 Omite los perfiles contaminados con señal de satelites,
2127 In: minHei = min_sat_range
2250 In: minHei = min_sat_range
@@ -2144,6 +2267,7 class CleanProfileSats(Operation):
2144 byRanges = False
2267 byRanges = False
2145 min_sats = None
2268 min_sats = None
2146 max_sats = None
2269 max_sats = None
2270 noise = 0
2147
2271
2148 def __init__(self, **kwargs):
2272 def __init__(self, **kwargs):
2149
2273
@@ -2177,16 +2301,19 class CleanProfileSats(Operation):
2177
2301
2178
2302
2179 def compareRanges(self,data, minHei,maxHei):
2303 def compareRanges(self,data, minHei,maxHei):
2180 ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref])
2304
2181 p_ref = 10*numpy.log10(ref.real)
2305 # ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref])
2182 m_ref = numpy.mean(p_ref)
2306 # p_ref = 10*numpy.log10(ref.real)
2307 # m_ref = numpy.mean(p_ref)
2308
2309 m_ref = self.noise
2183
2310
2184 sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
2311 sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
2185 p_sats = 10*numpy.log10(sats.real)
2312 p_sats = 10*numpy.log10(sats.real)
2186 m_sats = numpy.mean(p_sats)
2313 m_sats = numpy.mean(p_sats)
2187
2314
2188 if m_sats > (m_ref + self.th) and (m_sats > self.thdB):
2315 if m_sats > (m_ref + self.th): #and (m_sats > self.thdB):
2189 #print("msats: ",m_sats," mRef: ", m_ref, (m_sats - m_ref))
2316 #print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref))
2190 #print("Removing profiles...")
2317 #print("Removing profiles...")
2191 return False
2318 return False
2192
2319
@@ -2218,12 +2345,12 class CleanProfileSats(Operation):
2218 if not self.isConfig:
2345 if not self.isConfig:
2219 self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
2346 self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
2220 self.isConfig = True
2347 self.isConfig = True
2221
2348 #print(self.min_sats,self.max_sats)
2222 if dataOut.flagDataAsBlock:
2349 if dataOut.flagDataAsBlock:
2223 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
2350 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
2224
2351
2225 else:
2352 else:
2226
2353 self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref))
2227 if not self.isProfileClean(dataOut.data):
2354 if not self.isProfileClean(dataOut.data):
2228 return dataOut
2355 return dataOut
2229 #dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)
2356 #dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)
General Comments 0
You need to be logged in to leave comments. Login now