##// 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 67 double *sortdata = (double*)PyArray_DATA(data_array);
68 68 int lenOfData = (int)PyArray_SIZE(data_array) ;
69 double nums_min = lenOfData*0.75;
69 double nums_min = lenOfData*0.85;
70 70 if (nums_min <= 5) nums_min = 5;
71 71 double sump = 0;
72 72 double sumq = 0;
Show file before | Show file after | Hide comments
@@ -661,12 +661,10 class Project(Process):
661 661 elif ok == 'no_Read' and (not flag_no_read):
662 662 nProc_noRead = n_proc - 1
663 663 flag_no_read = True
664 print("No read")
665 664 continue
666 665 elif ok == 'new_Read':
667 666 nProc_noRead = 0
668 667 flag_no_read = False
669 print("read again")
670 668 continue
671 669 elif not ok:
672 670 #print("not ok",ok)
Show file before | Show file after | Hide comments
@@ -78,7 +78,7 def hildebrand_sekhon(data, navg):
78 78 sortdata = numpy.sort(data, axis=None)
79 79 '''
80 80 lenOfData = len(sortdata)
81 nums_min = lenOfData*0.2
81 nums_min = lenOfData*0.5
82 82
83 83 if nums_min <= 5:
84 84
@@ -106,6 +106,7 def hildebrand_sekhon(data, navg):
106 106 j += 1
107 107
108 108 lnoise = sump / j
109 return lnoise
109 110 '''
110 111 return _noise.hildebrand_sekhon(sortdata, navg)
111 112
@@ -380,7 +381,7 class Voltage(JROData):
380 381 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
381 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 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 392 if channel != None:
392 data = self.data[channel]
393 data = self.data[channel,ymin_index:ymax_index]
393 394 nChannels = 1
394 395 else:
395 data = self.data
396 data = self.data[:,ymin_index:ymax_index]
396 397 nChannels = self.nChannels
397 398
398 399 noise = numpy.zeros(nChannels)
@@ -407,10 +408,10 class Voltage(JROData):
407 408
408 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 413 if type == 1:
413 noise = self.getNoisebyHildebrand(channel)
414 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
414 415
415 416 return noise
416 417
@@ -437,6 +438,8 class Voltage(JROData):
437 438
438 439 class Spectra(JROData):
439 440
441 data_outlier = 0
442
440 443 def __init__(self):
441 444 '''
442 445 Constructor
@@ -475,6 +478,9 class Spectra(JROData):
475 478 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
476 479
477 480
481 self.max_nIncohInt = 1
482
483
478 484 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
479 485 """
480 486 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
@@ -482,12 +488,28 class Spectra(JROData):
482 488 Return:
483 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 508 noise = numpy.zeros(self.nChannels)
487 509 for channel in range(self.nChannels):
488 daux = self.data_spc[channel,
489 xmin_index:xmax_index, ymin_index:ymax_index]
490 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
510 daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index]
511
512 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt)
491 513
492 514 return noise
493 515
@@ -552,6 +574,7 class Spectra(JROData):
552 574 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
553 575 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
554 576
577
555 578 return normFactor
556 579
557 580 @property
@@ -582,7 +605,7 class Spectra(JROData):
582 605 def getPower(self):
583 606
584 607 factor = self.normFactor
585 z = self.data_spc / factor
608 z = numpy.divide(self.data_spc,factor)
586 609 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
587 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 58 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
59 59 data['spc'] = spc
60 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 62 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
63 63 if self.CODE == 'spc_moments':
64 64 data['moments'] = dataOut.moments
@@ -86,9 +86,10 class SpectraPlot(Plot):
86 86
87 87 data = self.data[-1]
88 88 z = data['spc']
89 print(z.shape, x.shape, y.shape)
89 #print(z.shape, x.shape, y.shape)
90 90 for n, ax in enumerate(self.axes):
91 noise = data['noise'][n]
91 #noise = data['noise'][n]
92 noise=62
92 93 if self.CODE == 'spc_moments':
93 94 mean = data['moments'][n, 1]
94 95 if ax.firsttime:
@@ -105,15 +106,15 class SpectraPlot(Plot):
105 106 if self.showprofile:
106 107 ax.plt_profile = self.pf_axes[n].plot(
107 108 data['rti'][n], y)[0]
108 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
109 color="k", linestyle="dashed", lw=1)[0]
109 # ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
110 # color="k", linestyle="dashed", lw=1)[0]
110 111 if self.CODE == 'spc_moments':
111 112 ax.plt_mean = ax.plot(mean, y, color='k')[0]
112 113 else:
113 114 ax.plt.set_array(z[n].T.ravel())
114 115 if self.showprofile:
115 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 118 if self.CODE == 'spc_moments':
118 119 ax.plt_mean.set_data(mean, y)
119 120 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
@@ -274,7 +275,7 class RTIPlot(Plot):
274 275
275 276 self.x = self.data.times
276 277 self.y = self.data.yrange
277
278 #print(" x, y: ",self.x, self.y)
278 279 self.z = self.data[self.CODE]
279 280 self.z = numpy.array(self.z, dtype=float)
280 281 self.z = numpy.ma.masked_invalid(self.z)
@@ -816,6 +817,7 class NoiselessSpectraPlot(Plot):
816 817 plot_type = 'pcolor'
817 818 buffering = False
818 819 channelList = []
820 last_noise = None
819 821
820 822 def setup(self):
821 823
@@ -842,18 +844,22 class NoiselessSpectraPlot(Plot):
842 844 self.update_list(dataOut)
843 845 data = {}
844 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
854 data['rti'] = dataOut.getPower() - n1
851 if self.last_noise == None:
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 863 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
858 864
859 865 return data, meta
@@ -877,7 +883,7 class NoiselessSpectraPlot(Plot):
877 883 z = data['spc']
878 884
879 885 for n, ax in enumerate(self.axes):
880 noise = data['noise'][n]
886 #noise = data['noise'][n]
881 887
882 888 if ax.firsttime:
883 889 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
@@ -893,16 +899,15 class NoiselessSpectraPlot(Plot):
893 899 if self.showprofile:
894 900 ax.plt_profile = self.pf_axes[n].plot(
895 901 data['rti'][n], y)[0]
896 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
897 color="k", linestyle="dashed", lw=1)[0]
902
898 903
899 904 else:
900 905 ax.plt.set_array(z[n].T.ravel())
901 906 if self.showprofile:
902 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 913 class NoiselessRTIPlot(Plot):
@@ -917,6 +922,7 class NoiselessRTIPlot(Plot):
917 922 channelList = []
918 923 elevationList = []
919 924 azimuthList = []
925 last_noise = None
920 926
921 927 def setup(self):
922 928 self.xaxis = 'time'
@@ -945,19 +951,21 class NoiselessRTIPlot(Plot):
945 951 data = {}
946 952 meta = {}
947 953
948 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
949 (nch, nff, nh) = dataOut.data_spc.shape
950 #print(nch, nff, nh)
951 if nch != 1:
952 aux = []
953 for c in self.channelList:
954 aux.append(n0[c])
955 n0 = numpy.asarray(aux)
956 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
957 #print(dataOut.elevationList, dataOut.azimuthList)
958 #print(dataOut.channelList)
959 data['noiseless_rti'] = dataOut.getPower() - noise
960 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
954
955 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
956 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
957 #print("noise: ",n0, dataOut.normFactor, norm, dataOut.nIncohInt, dataOut.max_nIncohInt)
958 if self.last_noise == None:
959 self.last_noise = n0
960 else:
961 n0 = (n0*0.2 + self.last_noise*0.8)
962 self.last_noise = n0
963
964
965 data['noiseless_rti'] = dataOut.getPower() - n0
966
967 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
968 #print(noise)
961 969 return data, meta
962 970
963 971 def plot(self):
@@ -986,6 +994,7 class NoiselessRTIPlot(Plot):
986 994 else:
987 995 x, y, z = self.fill_gaps(*self.decimate())
988 996 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
997 #print("plot shapes ", z.shape, x.shape, y.shape)
989 998 for n, ax in enumerate(self.axes):
990 999 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
991 1000 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
@@ -999,9 +1008,6 class NoiselessRTIPlot(Plot):
999 1008 if self.showprofile:
1000 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 1011 else:
1006 1012 ax.collections.remove(ax.collections[0])
1007 1013 ax.plt = ax.pcolormesh(x, y, z[n].T,
@@ -1011,6 +1017,164 class NoiselessRTIPlot(Plot):
1011 1017 )
1012 1018 if self.showprofile:
1013 1019 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1014 if "noise" in self.data:
1015 ax.plot_noise.set_data(numpy.repeat(
1016 data['noise'][n], len(self.y)), self.y)
1020 # if "noise" in self.data:
1021 # #ax.plot_noise.set_data(numpy.repeat(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 185 self.nsa = self.nsamplesPulse[0,0] #ngates
186 186 self.nchannels = len(self.beamCode)
187 187 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
188 #print("IPPS secs: ",self.ippSeconds)
188 189 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
189 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 20 from schainpy.model.data import _noise
21 21
22 22 from schainpy.utils import log
23
23 import matplotlib.pyplot as plt
24 24 from scipy.optimize import curve_fit
25 25
26 26 class SpectraProc(ProcessingUnit):
@@ -125,7 +125,7 class SpectraProc(ProcessingUnit):
125 125 self.dataOut.flagShiftFFT = False
126 126
127 127 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
128
128 #print("run spc proc")
129 129 if self.dataIn.type == "Spectra":
130 130
131 131 try:
@@ -199,6 +199,7 class SpectraProc(ProcessingUnit):
199 199 if self.profIndex == nProfiles:
200 200
201 201 self.__updateSpecFromVoltage()
202
202 203 if pairsList == None:
203 204 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
204 205 else:
@@ -208,6 +209,7 class SpectraProc(ProcessingUnit):
208 209 self.firstdatatime = None
209 210 self.profIndex = 0
210 211
212
211 213 else:
212 214 raise ValueError("The type of input object '%s' is not valid".format(
213 215 self.dataIn.type))
@@ -529,55 +531,56 class getNoise(Operation):
529 531 # validacion de velocidades
530 532 indminPoint = None
531 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:
538 minFreq = freqrange[0]
542 if maxFreq == None:
543 maxFreq = freqrange[-1]
539 544
540 if maxFreq == None:
541 maxFreq = freqrange[-1]
545 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
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):
544 if self.warnings:
545 print('minFreq: %.2f is out of the frequency range' % (minFreq))
546 print('minFreq is setting to %.2f' % (freqrange[0]))
547 minFreq = freqrange[0]
551 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
552 if self.warnings:
553 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
554 print('maxFreq is setting to %.2f' % (freqrange[-1]))
555 maxFreq = freqrange[-1]
548 556
549 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
550 if self.warnings:
551 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
552 print('maxFreq is setting to %.2f' % (freqrange[-1]))
553 maxFreq = freqrange[-1]
557 indminPoint = numpy.where(freqrange >= minFreq)
558 indmaxPoint = numpy.where(freqrange <= maxFreq)
554 559
555 indminPoint = numpy.where(freqrange >= minFreq)
556 indmaxPoint = numpy.where(freqrange <= maxFreq)
560 else:
557 561
558 else:
559 velrange = self.dataOut.getVelRange(1)
562 velrange = self.dataOut.getVelRange(1)
560 563
561 if minVel == None:
562 minVel = velrange[0]
564 if minVel == None:
565 minVel = velrange[0]
563 566
564 if maxVel == None:
565 maxVel = velrange[-1]
567 if maxVel == None:
568 maxVel = velrange[-1]
566 569
567 if (minVel < velrange[0]) or (minVel > maxVel):
568 if self.warnings:
569 print('minVel: %.2f is out of the velocity range' % (minVel))
570 print('minVel is setting to %.2f' % (velrange[0]))
571 minVel = velrange[0]
570 if (minVel < velrange[0]) or (minVel > maxVel):
571 if self.warnings:
572 print('minVel: %.2f is out of the velocity range' % (minVel))
573 print('minVel is setting to %.2f' % (velrange[0]))
574 minVel = velrange[0]
572 575
573 if (maxVel > velrange[-1]) or (maxVel < minVel):
574 if self.warnings:
575 print('maxVel: %.2f is out of the velocity range' % (maxVel))
576 print('maxVel is setting to %.2f' % (velrange[-1]))
577 maxVel = velrange[-1]
576 if (maxVel > velrange[-1]) or (maxVel < minVel):
577 if self.warnings:
578 print('maxVel: %.2f is out of the velocity range' % (maxVel))
579 print('maxVel is setting to %.2f' % (velrange[-1]))
580 maxVel = velrange[-1]
578 581
579 indminPoint = numpy.where(velrange >= minVel)
580 indmaxPoint = numpy.where(velrange <= maxVel)
582 indminPoint = numpy.where(velrange >= minVel)
583 indmaxPoint = numpy.where(velrange <= maxVel)
581 584
582 585
583 586 # seleccion de indices para rango
@@ -606,23 +609,30 class getNoise(Operation):
606 609 maxIndex = self.dataOut.nHeights - 1
607 610 #############################################################3
608 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:
611 minIndexFFT = indminPoint[0][0]
612 except:
613 minIndexFFT = 0
618 try:
619 maxIndexFFT = indmaxPoint[0][-1]
620 except:
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 624 self.dataOut.noise_estimation = None
622 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
623
625 noise = None
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 633 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
625 634 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
635
626 636 return self.dataOut
627 637
628 638
@@ -705,7 +715,7 class CleanRayleigh(Operation):
705 715
706 716
707 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 719 if not self.isConfig :
710 720
711 721 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
@@ -736,7 +746,10 class CleanRayleigh(Operation):
736 746 if numpy.any(jspc) :
737 747 #print( jspc.shape, jcspc.shape)
738 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 753 self.__dataReady = False
741 754 #print( jspc.shape, jcspc.shape)
742 755 dataOut.flagNoData = False
@@ -748,8 +761,11 class CleanRayleigh(Operation):
748 761 #print( len(self.buffer))
749 762 if numpy.any(self.buffer):
750 763 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
751 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
752 self.buffer3 += dataOut.data_dc
764 try:
765 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
766 self.buffer3 += dataOut.data_dc
767 except:
768 pass
753 769 else:
754 770 self.buffer = dataOut.data_spc
755 771 self.buffer2 = dataOut.data_cspc
@@ -762,7 +778,9 class CleanRayleigh(Operation):
762 778
763 779
764 780 #index = tini.tm_hour*12+tini.tm_min/5
765 '''REVISAR'''
781 '''
782 REVISAR
783 '''
766 784 # jspc = jspc/self.nFFTPoints/self.normFactor
767 785 # jcspc = jcspc/self.nFFTPoints/self.normFactor
768 786
@@ -776,6 +794,7 class CleanRayleigh(Operation):
776 794
777 795 dataOut.data_dc = self.buffer3
778 796 dataOut.nIncohInt *= self.nIntProfiles
797 dataOut.max_nIncohInt = self.nIntProfiles
779 798 dataOut.utctime = self.currentTime #tiempo promediado
780 799 #print("Time: ",time.localtime(dataOut.utctime))
781 800 # dataOut.data_spc = sat_spectra
@@ -787,7 +806,7 class CleanRayleigh(Operation):
787 806 return dataOut
788 807
789 808 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
790 #print("OP cleanRayleigh")
809 print("OP cleanRayleigh")
791 810 #import matplotlib.pyplot as plt
792 811 #for k in range(149):
793 812 #channelsProcssd = []
@@ -803,6 +822,8 class CleanRayleigh(Operation):
803 822
804 823 ###ONLY FOR TEST:
805 824 raxs = math.ceil(math.sqrt(self.nPairs))
825 if raxs == 0:
826 raxs = 1
806 827 caxs = math.ceil(self.nPairs/raxs)
807 828 if self.nPairs <4:
808 829 raxs = 2
@@ -1100,12 +1121,13 class IntegrationFaradaySpectra(Operation):
1100 1121 __dataReady = False
1101 1122
1102 1123 __timeInterval = None
1103
1124 n_ints = None #matriz de numero de integracions (CH,HEI)
1104 1125 n = None
1105 1126 minHei_ind = None
1106 1127 maxHei_ind = None
1107 1128 navg = 1.0
1108 1129 factor = 0.0
1130 dataoutliers = None # (CHANNELS, HEIGHTS)
1109 1131
1110 1132 def __init__(self):
1111 1133
@@ -1138,6 +1160,7 class IntegrationFaradaySpectra(Operation):
1138 1160 self.navg = avg
1139 1161 #self.ByLags = dataOut.ByLags ###REDEFINIR
1140 1162 self.ByLags = False
1163 self.maxProfilesInt = 1
1141 1164
1142 1165 if DPL != None:
1143 1166 self.DPL=DPL
@@ -1251,6 +1274,8 class IntegrationFaradaySpectra(Operation):
1251 1274 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1252 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 1279 for k in range(self.minHei_ind,self.maxHei_ind):
1255 1280 try:
1256 1281 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
@@ -1266,22 +1291,31 class IntegrationFaradaySpectra(Operation):
1266 1291 #sortIDs=[]
1267 1292 outliers_IDs=[]
1268 1293
1269 for j in range(self.nProfiles):
1294 for j in range(self.nProfiles): #frecuencias en el tiempo
1270 1295 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1271 1296 # continue
1272 1297 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1273 1298 # continue
1274 1299 buffer=buffer1[:,j]
1275 1300 sortdata = numpy.sort(buffer, axis=None)
1301
1276 1302 sortID=buffer.argsort()
1277 1303 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1278 1304
1279 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 1314 indexes.append(index)
1282 1315 #sortIDs.append(sortID)
1283 1316 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1284 1317
1318 #print("Outliers: ",outliers_IDs)
1285 1319 outliers_IDs=numpy.array(outliers_IDs)
1286 1320 outliers_IDs=outliers_IDs.ravel()
1287 1321 outliers_IDs=numpy.unique(outliers_IDs)
@@ -1289,31 +1323,45 class IntegrationFaradaySpectra(Operation):
1289 1323 indexes=numpy.array(indexes)
1290 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 1337 if indexmin != buffer1.shape[0]:
1293 1338 if self.nChannels > 1:
1294 1339 cspc_outliers_exist= True
1295 1340
1296 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 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 1350 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1303 ###cspc IDs
1304 #indexmin_cspc+=indexmin_cspc
1351
1305 1352 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1306 1353
1307 #if not breakFlag:
1354
1355
1308 1356 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1309 1357 if cspc_outliers_exist :
1310 #sortdata=numpy.sort(buffer_cspc,axis=0)
1311 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1358
1312 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 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 1366 try:
1319 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 1378 buffer=None
1330 1379 bufferH=None
1331 1380 buffer1=None
1332 1381 buffer_cspc=None
1333 1382
1334 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1335 #exit()
1336 1383
1337 1384 buffer=None
1338 #print(self.__buffer_spc[:,1,3,20,0])
1339 #print(self.__buffer_spc[:,1,5,37,0])
1340 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1385
1386 #data_spc = numpy.sum(self.__buffer_spc,axis=0)
1387 data_spc = numpy.nansum(self.__buffer_spc,axis=0)
1341 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 1391 except:
1344 1392 #print("No cpsc",e)
1345 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 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 1401 self.__buffer_spc = []
1357 1402 self.__buffer_cspc = []
1358 1403 self.__buffer_dc = 0
1359 1404 self.__profIndex = 0
1360 #print("pushData Done")
1405
1361 1406 return data_spc, data_cspc, data_dc, n
1362 1407
1363 1408 def byProfiles(self, *args):
@@ -1372,7 +1417,7 class IntegrationFaradaySpectra(Operation):
1372 1417 if self.__profIndex == self.n:
1373 1418
1374 1419 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1375 self.n = n
1420 self.n_ints = n
1376 1421 self.__dataReady = True
1377 1422
1378 1423 return avgdata_spc, avgdata_cspc, avgdata_dc
@@ -1388,7 +1433,7 class IntegrationFaradaySpectra(Operation):
1388 1433
1389 1434 if (datatime - self.__initime) >= self.__integrationtime:
1390 1435 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1391 self.n = n
1436 self.n_ints = n
1392 1437 self.__dataReady = True
1393 1438
1394 1439 return avgdata_spc, avgdata_cspc, avgdata_dc
@@ -1450,7 +1495,7 class IntegrationFaradaySpectra(Operation):
1450 1495 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
1451 1496 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1452 1497 self.dataOut.data_dc = avgdata_dc
1453
1498 self.dataOut.data_outlier = self.dataOutliers
1454 1499
1455 1500 else:
1456 1501 self.dataOut.dataLag_spc = avgdata_spc
@@ -1462,10 +1507,12 class IntegrationFaradaySpectra(Operation):
1462 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 1513 self.dataOut.utctime = avgdatatime
1467 1514 self.dataOut.flagNoData = False
1468
1515 #print("Faraday Integration DONE...")
1469 1516 return self.dataOut
1470 1517
1471 1518 class removeInterference(Operation):
@@ -1705,7 +1752,8 class IncohInt(Operation):
1705 1752 __dataReady = False
1706 1753
1707 1754 __timeInterval = None
1708
1755 incohInt = 0
1756 nOutliers = 0
1709 1757 n = None
1710 1758
1711 1759 def __init__(self):
@@ -1734,7 +1782,8 class IncohInt(Operation):
1734 1782 self.__profIndex = 0
1735 1783 self.__dataReady = False
1736 1784 self.__byTime = False
1737
1785 self.incohInt = 0
1786 self.nOutliers = 0
1738 1787 if n is None and timeInterval is None:
1739 1788 raise ValueError("n or timeInterval should be specified ...")
1740 1789
@@ -1788,7 +1837,7 class IncohInt(Operation):
1788 1837 self.__buffer_spc = 0
1789 1838 self.__buffer_cspc = 0
1790 1839 self.__buffer_dc = 0
1791 self.__profIndex = 0
1840
1792 1841
1793 1842 return data_spc, data_cspc, data_dc, n
1794 1843
@@ -1845,6 +1894,9 class IncohInt(Operation):
1845 1894 if n == 1:
1846 1895 return dataOut
1847 1896
1897 if dataOut.flagNoData == True:
1898 return dataOut
1899
1848 1900 dataOut.flagNoData = True
1849 1901
1850 1902 if not self.isConfig:
@@ -1855,15 +1907,21 class IncohInt(Operation):
1855 1907 dataOut.data_spc,
1856 1908 dataOut.data_cspc,
1857 1909 dataOut.data_dc)
1858
1910 self.incohInt += dataOut.nIncohInt
1911 self.nOutliers += dataOut.data_outlier
1859 1912 if self.__dataReady:
1860 1913
1861 1914 dataOut.data_spc = avgdata_spc
1862 1915 dataOut.data_cspc = avgdata_cspc
1863 1916 dataOut.data_dc = avgdata_dc
1864 dataOut.nIncohInt *= self.n
1917 dataOut.nIncohInt = self.incohInt
1918 dataOut.data_outlier = self.nOutliers
1865 1919 dataOut.utctime = avgdatatime
1866 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 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 8 from time import time
9 9 import datetime
10 10 import numpy
11 import copy
11 #import copy
12 from schainpy.model.data import _noise
12 13
13 14 class VoltageProc(ProcessingUnit):
14 15
@@ -1759,6 +1760,8 class SSheightProfiles2(Operation):
1759 1760 profileShape = None
1760 1761 sshProfiles = None
1761 1762 profileIndex = None
1763 deltaHeight = None
1764 init_range = None
1762 1765
1763 1766 def __init__(self, **kwargs):
1764 1767
@@ -1780,7 +1783,8 class SSheightProfiles2(Operation):
1780 1783 if residue != 0:
1781 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 1788 #numberProfile = self.nsamples
1785 1789 self.numberSamples = (self.__nHeis - self.nsamples)/self.step
1786 1790
@@ -1800,21 +1804,25 class SSheightProfiles2(Operation):
1800 1804
1801 1805 if repeat is not None:
1802 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 1810 for i in range(int(self.new_nHeights)): #nuevas alturas
1805 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 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 1816 for j in range(self.__nChannels): #en los cananles
1811 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1812
1817 self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
1818 #print("new profs Done")
1813 1819
1814 1820
1815 1821
1816 1822 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1817 1823
1824 if dataOut.flagNoData == True:
1825 return dataOut
1818 1826 dataOut.flagNoData = True
1819 1827 #print("init data shape:", dataOut.data.shape)
1820 1828 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
@@ -1838,24 +1846,26 class SSheightProfiles2(Operation):
1838 1846 #print("dataOut nProfiles:", dataOut.nProfiles)
1839 1847 for profile in range(nprof):
1840 1848 if dataOut.flagDataAsBlock:
1849 #print("read blocks")
1841 1850 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
1842 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 1854 if profile == 0:
1845 1855 dataBlock = self.sshProfiles.copy()
1846 1856 else: #by blocks
1847 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 1860 profileIndex = self.nsamples
1851 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1852 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1861 #deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1862 ippSeconds = (self.deltaHeight*1.0e-6)/(0.15)
1853 1863
1854 1864
1855 1865 dataOut.data = dataBlock
1856 dataOut.heightList = numpy.arange(self.new_nHeights) *self.step*deltaHeight + dataOut.heightList[0]
1857 #print("show me: ",dataOut.nHeights, self.new_nHeights)
1858 #dataOut.nHeights = int(self.new_nHeights)
1866 #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
1867 dataOut.heightList = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
1868
1859 1869 dataOut.ippSeconds = ippSeconds
1860 1870 dataOut.step = self.step
1861 1871 dataOut.flagNoData = False
@@ -1877,17 +1887,17 class SSheightProfiles2(Operation):
1877 1887
1878 1888 #import skimage.color
1879 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 1896 isConfig = False
1887 1897 n = None
1888 __timeInterval = None
1898
1889 1899 __profIndex = 0
1890 __byTime = False
1900
1891 1901 __dataReady = False
1892 1902 __buffer_data = []
1893 1903 __buffer_times = []
@@ -1901,6 +1911,7 class clean2DArray(Operation):
1901 1911 end_prof = 0
1902 1912 n_profiles = 0
1903 1913 first_utcBlock = None
1914 outliers_IDs_list = []
1904 1915 __dh = 0
1905 1916
1906 1917 def __init__(self, **kwargs):
@@ -1908,42 +1919,122 class clean2DArray(Operation):
1908 1919 Operation.__init__(self, **kwargs)
1909 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 1924 if n == None and timeInterval == None:
1914 1925 raise ValueError("nprofiles or timeInterval should be specified ...")
1915 1926
1916 1927 if n != None:
1917 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 1933 self.__profIndex = 0
1925 1934 self.buffer = None
1926 1935 self.lenProfileOut = 1
1927
1936 self.n_prof_released = 0
1937 self.heightList = dataOut.heightList
1928 1938 self.init_prof = 0
1929 1939 self.end_prof = 0
1930 1940 self.n_profiles = 0
1931 1941 self.first_utcBlock = None
1932 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 2025 def cleanOutliersByBlock(self):
1935 2026 #print(self.__buffer_data[0].shape)
1936 2027 data = self.__buffer_data#.copy()
1937 2028 #print("cleaning shape inpt: ",data.shape)
1938
2029 '''
1939 2030 self.__buffer_data = []
1940 2031
1941 2032 spectrum = numpy.fft.fft2(data, axes=(0,2))
1942 2033 #print("spc : ",spectrum.shape)
1943 2034 (nch,nsamples, nh) = spectrum.shape
1944
1945 #nch =
1946
2035 data2 = None
2036 #print(data.shape)
2037 spectrum2 = spectrum.copy()
1947 2038 for ch in range(nch):
1948 2039 dh = self.__dh
1949 2040 dt1 = (dh*1.0e-6)/(0.15)
@@ -1951,57 +2042,48 class clean2DArray(Operation):
1951 2042 freqv = numpy.fft.fftfreq(nh, d=dt1)
1952 2043 freqh = numpy.fft.fftfreq(self.n, d=dt2)
1953 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 2080 #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
2081 '''
2082 #print("cleanOutliersByBlock Done")
1976 2083
1977 print("cleanOutliersByBlock Done")
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
2084 return self.filterSatsProfiles()
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 2088 def fillBuffer(self, data, datatime):
2007 2089
@@ -2011,54 +2093,60 class clean2DArray(Operation):
2011 2093 else:
2012 2094 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2013 2095 self.__profIndex += 1
2014 self.__buffer_times.append(datatime)
2096 #self.__buffer_times.append(datatime)
2015 2097
2016 2098 def getData(self, data, datatime=None):
2017 2099
2018 2100 if self.__profIndex == 0:
2019 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 2117 if dataBlock is None:
2029 return None
2118 return None, None
2030 2119
2031 2120
2032 2121
2033 2122 return dataBlock
2034 2123
2035 def releaseBlock(self, dataOut):
2124 def releaseBlock(self):
2036 2125
2037 2126 if self.n % self.lenProfileOut != 0:
2038 2127 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2039 2128 return None
2040 2129
2041 dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2042 dataOut.profileIndex = self.lenProfileOut
2043 dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds
2130 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
2131
2044 2132 self.init_prof = self.end_prof
2045 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 2142 #print("run op buffer 2D")
2057 2143 self.nChannels = dataOut.nChannels
2058 2144 self.nHeights = dataOut.nHeights
2059 2145
2060 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 2150 self.isConfig = True
2063 2151
2064 2152 dataBlock = None
@@ -2066,6 +2154,7 class clean2DArray(Operation):
2066 2154 if not dataOut.buffer_empty: #hay datos acumulados
2067 2155
2068 2156 if self.init_prof == 0:
2157 self.n_prof_released = 0
2069 2158 self.lenProfileOut = nProfilesOut
2070 2159 dataOut.flagNoData = False
2071 2160 #print("tp 2 ",dataOut.data.shape)
@@ -2074,13 +2163,45 class clean2DArray(Operation):
2074 2163 #print("rel: ",self.buffer[:,-1,:])
2075 2164 self.init_prof = 0
2076 2165 self.end_prof = self.lenProfileOut
2077 self.first_utcBlock = dataOut.utctime
2166
2078 2167 dataOut.nProfiles = self.lenProfileOut
2079 dataOut.error = False
2168 if nProfilesOut == 1:
2169 dataOut.flagDataAsBlock = False
2170 else:
2080 2171 dataOut.flagDataAsBlock = True
2081 2172 #print("prof: ",self.init_prof)
2082 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 2207 #print("tp 223 ",dataOut.data.shape)
@@ -2098,13 +2219,16 class clean2DArray(Operation):
2098 2219 sys.exit()
2099 2220
2100 2221 if self.__dataReady:
2222 #print("omitting: ", len(self.outliers_IDs_list))
2101 2223 self.__count_exec = 0
2102 2224 #dataOut.data =
2103 self.buffer = numpy.flip(dataBlock, axis=1)
2104
2225 #self.buffer = numpy.flip(dataBlock, axis=1)
2226 self.buffer = dataBlock
2227 self.first_utcBlock = self.__initime
2105 2228 dataOut.utctime = self.__initime
2106 2229 dataOut.nProfiles = self.__profIndex
2107 2230 #dataOut.flagNoData = False
2231 self.init_prof = 0
2108 2232 self.__profIndex = 0
2109 2233 self.__initime = None
2110 2234 dataBlock = None
@@ -2112,16 +2236,15 class clean2DArray(Operation):
2112 2236 dataOut.error = False
2113 2237 dataOut.useInputBuffer = True
2114 2238 dataOut.buffer_empty = False
2115 print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
2116 int(dataOut.nProfiles),int(dataOut.nHeights)))
2117 #return None
2239 #print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
2240
2118 2241
2119 2242
2120 2243 #print(self.__count_exec)
2121 2244
2122 2245 return dataOut
2123 2246
2124 class CleanProfileSats(Operation):
2247 class RemoveProfileSats(Operation):
2125 2248 '''
2126 2249 Omite los perfiles contaminados con señal de satelites,
2127 2250 In: minHei = min_sat_range
@@ -2144,6 +2267,7 class CleanProfileSats(Operation):
2144 2267 byRanges = False
2145 2268 min_sats = None
2146 2269 max_sats = None
2270 noise = 0
2147 2271
2148 2272 def __init__(self, **kwargs):
2149 2273
@@ -2177,16 +2301,19 class CleanProfileSats(Operation):
2177 2301
2178 2302
2179 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])
2181 p_ref = 10*numpy.log10(ref.real)
2182 m_ref = numpy.mean(p_ref)
2304
2305 # ref = data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_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 2311 sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
2185 2312 p_sats = 10*numpy.log10(sats.real)
2186 2313 m_sats = numpy.mean(p_sats)
2187 2314
2188 if m_sats > (m_ref + self.th) and (m_sats > self.thdB):
2189 #print("msats: ",m_sats," mRef: ", m_ref, (m_sats - m_ref))
2315 if m_sats > (m_ref + self.th): #and (m_sats > self.thdB):
2316 #print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref))
2190 2317 #print("Removing profiles...")
2191 2318 return False
2192 2319
@@ -2218,12 +2345,12 class CleanProfileSats(Operation):
2218 2345 if not self.isConfig:
2219 2346 self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
2220 2347 self.isConfig = True
2221
2348 #print(self.min_sats,self.max_sats)
2222 2349 if dataOut.flagDataAsBlock:
2223 2350 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
2224 2351
2225 2352 else:
2226
2353 self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref))
2227 2354 if not self.isProfileClean(dataOut.data):
2228 2355 return dataOut
2229 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