##// END OF EJS Templates
schain
RSS
grafico outliers, grafico de integraciones, funciones de obtención ruído, etc
joabAM -
r1540:c8092aa564c5
parent child
Show More
Side by Side Unified
Context file:
r1540:c8092aa564c5 Loading diff...:
Show file before | Show file after | Hide comments
@@ -20,7 +20,7 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
20 20 }
21 21 double *sortdata = (double*)PyArray_DATA(data_array);
22 22 int lenOfData = (int)PyArray_SIZE(data_array) ;
23 double nums_min = lenOfData*0.2;
23 double nums_min = lenOfData*0.5;//0.2
24 24 if (nums_min <= 5) nums_min = 5;
25 25 double sump = 0;
26 26 double sumq = 0;
Show file before | Show file after | Hide comments
@@ -340,7 +340,7 class Project(Process):
340 340 idList = list(self.configurations.keys())
341 341 id = int(self.id) * 10
342 342
343 while True:
343 while 1:
344 344 id += 1
345 345
346 346 if str(id) in idList:
@@ -633,7 +633,7 class Project(Process):
633 633 n = len(self.configurations)
634 634 flag_no_read = False
635 635 nProc_noRead = 0
636 while not err:
636 while 1:
637 637 #print("STAR")
638 638 n_proc = 0
639 639 for conf in self.getUnits():
@@ -672,6 +672,7 class Project(Process):
672 672
673 673 if n == 0:
674 674 err = True
675 break
675 676
676 677 def run(self):
677 678
Show file before | Show file after | Hide comments
@@ -199,6 +199,7 class JROData(GenericData):
199 199 codeList = []
200 200 azimuthList = []
201 201 elevationList = []
202 last_noise = None
202 203
203 204 def __str__(self):
204 205
@@ -641,7 +642,7 class Spectra(JROData):
641 642
642 643 def setValue(self, value):
643 644
644 print("This property should not be initialized")
645 print("This property should not be initialized", value)
645 646
646 647 return
647 648
@@ -875,6 +876,7 class Parameters(Spectra):
875 876 data_param = None # Parameters obtained
876 877 data_pre = None # Data Pre Parametrization
877 878 data_SNR = None # Signal to Noise Ratio
879 data_outlier = None
878 880 abscissaList = None # Abscissa, can be velocities, lags or time
879 881 utctimeInit = None # Initial UTC time
880 882 paramInterval = None # Time interval to calculate Parameters in seconds
@@ -889,7 +891,7 class Parameters(Spectra):
889 891 nAvg = None
890 892 noise_estimation = None
891 893 GauSPC = None # Fit gaussian SPC
892
894 max_nIncohInt = 1
893 895 def __init__(self):
894 896 '''
895 897 Constructor
Show file before | Show file after | Hide comments
@@ -580,7 +580,7 class Plot(Operation):
580 580
581 581 self.sender_queue.append(last_time)
582 582
583 while True:
583 while 1:
584 584 try:
585 585 tm = self.sender_queue.popleft()
586 586 except IndexError:
Show file before | Show file after | Hide comments
@@ -53,7 +53,7 class SnrPlot(RTIPlot):
53 53
54 54 meta = {}
55 55 data = {
56 'snr': 10 * numpy.log10(dataOut.data_snr)
56 'snr': 10 * numpy.log10(dataOut.data_snr * dataOut.nIncohInt/ dataOut.max_nIncohInt)
57 57 }
58 58 #print(data['snr'])
59 59 return data, meta
Show file before | Show file after | Hide comments
@@ -58,7 +58,10 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 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
62 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
63 data['noise'] = noise[0]
64
62 65 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
63 66 if self.CODE == 'spc_moments':
64 67 data['moments'] = dataOut.moments
@@ -88,8 +91,8 class SpectraPlot(Plot):
88 91 z = data['spc']
89 92 #print(z.shape, x.shape, y.shape)
90 93 for n, ax in enumerate(self.axes):
91 #noise = data['noise'][n]
92 noise=62
94 noise = self.data['noise'][n]
95 #print(noise)
93 96 if self.CODE == 'spc_moments':
94 97 mean = data['moments'][n, 1]
95 98 if ax.firsttime:
@@ -268,7 +271,11 class RTIPlot(Plot):
268 271 data = {}
269 272 meta = {}
270 273 data['rti'] = dataOut.getPower()
271 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
274
275 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
276 noise = 10*numpy.log10(dataOut.getNoise()/float(norm))
277 data['noise'] = noise
278
272 279 return data, meta
273 280
274 281 def plot(self):
@@ -327,9 +334,7 class RTIPlot(Plot):
327 334 if self.showprofile:
328 335 ax.plot_profile.set_data(data[self.CODE][n], self.y)
329 336 if "noise" in self.data:
330
331 ax.plot_noise.set_data(numpy.repeat(
332 data['noise'][n], len(self.y)), self.y)
337 ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
333 338
334 339
335 340 class CoherencePlot(RTIPlot):
@@ -401,12 +406,17 class NoisePlot(Plot):
401 406 self.titles = ['Noise']
402 407 self.colorbar = False
403 408 self.plots_adjust.update({'right': 0.85 })
409 #if not self.titles:
410 self.titles = ['Noise Plot']
404 411
405 412 def update(self, dataOut):
406 413
407 414 data = {}
408 415 meta = {}
409 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
416 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
417 #noise = 10*numpy.log10(dataOut.getNoise()/norm)
418 noise = 10*numpy.log10(dataOut.getNoise())
419 noise = noise.reshape(dataOut.nChannels, 1)
410 420 data['noise'] = noise
411 421 meta['yrange'] = numpy.array([])
412 422
@@ -491,11 +501,11 class SpectraCutPlot(Plot):
491 501 self.nplots = len(self.data.channels)
492 502 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
493 503 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
494 self.width = 4.2 * self.ncols + 2.5
504 self.width = 4.5 * self.ncols + 2.5
495 505 self.height = 4.8 * self.nrows
496 506 self.ylabel = 'Power [dB]'
497 507 self.colorbar = False
498 self.plots_adjust.update({'left':0.05, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
508 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
499 509
500 510 if len(self.selectedHeightsList) > 0:
501 511 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
@@ -519,7 +529,11 class SpectraCutPlot(Plot):
519 529 #print(self.height_index)
520 530 data = {}
521 531 meta = {}
522 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
532
533 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
534 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
535
536 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - n0
523 537
524 538 data['spc'] = spc
525 539 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
@@ -824,14 +838,15 class NoiselessSpectraPlot(Plot):
824 838 self.nplots = len(self.data.channels)
825 839 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
826 840 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
827 self.height = 2.6 * self.nrows
841 self.height = 3.5 * self.nrows
828 842
829 843 self.cb_label = 'dB'
830 844 if self.showprofile:
831 self.width = 4 * self.ncols
845 self.width = 5.8 * self.ncols
832 846 else:
833 self.width = 3.5 * self.ncols
834 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
847 self.width = 4.8* self.ncols
848 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
849
835 850 self.ylabel = 'Range [km]'
836 851
837 852
@@ -848,11 +863,6 class NoiselessSpectraPlot(Plot):
848 863 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
849 864 n0 = 10*numpy.log10(dataOut.getNoise()/float(norm))
850 865
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 866
857 867 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
858 868
@@ -954,18 +964,11 class NoiselessRTIPlot(Plot):
954 964
955 965 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
956 966 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 967
968 data['noise'] = n0[0]
964 969
965 970 data['noiseless_rti'] = dataOut.getPower() - n0
966 971
967 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
968 #print(noise)
969 972 return data, meta
970 973
971 974 def plot(self):
@@ -976,6 +979,7 class NoiselessRTIPlot(Plot):
976 979 self.z = numpy.array(self.z, dtype=float)
977 980 self.z = numpy.ma.masked_invalid(self.z)
978 981
982
979 983 try:
980 984 if self.channelList != None:
981 985 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
@@ -989,6 +993,8 class NoiselessRTIPlot(Plot):
989 993
990 994 self.titles = ['{} Channel {}'.format(
991 995 self.CODE.upper(), x) for x in self.channelList]
996
997
992 998 if self.decimation is None:
993 999 x, y, z = self.fill_gaps(self.x, self.y, self.z)
994 1000 else:
@@ -996,6 +1002,8 class NoiselessRTIPlot(Plot):
996 1002 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
997 1003 #print("plot shapes ", z.shape, x.shape, y.shape)
998 1004 for n, ax in enumerate(self.axes):
1005
1006
999 1007 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1000 1008 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1001 1009 data = self.data[-1]
@@ -1027,14 +1035,14 class OutliersRTIPlot(Plot):
1027 1035 Plot for data_xxxx object
1028 1036 '''
1029 1037
1030 CODE = 'outlier'
1038 CODE = 'outlier_rtc' # Range Time Counts
1031 1039 colormap = 'cool'
1032 1040 plot_type = 'pcolorbuffer'
1033 1041
1034 1042 def setup(self):
1035 1043 self.xaxis = 'time'
1036 1044 self.ncols = 1
1037 self.nrows = self.data.shape('outlier')[0]
1045 self.nrows = self.data.shape('outlier_rtc')[0]
1038 1046 self.nplots = self.nrows
1039 1047 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1040 1048
@@ -1049,7 +1057,7 class OutliersRTIPlot(Plot):
1049 1057 def update(self, dataOut):
1050 1058
1051 1059 data = {}
1052 data['outlier'] = dataOut.data_outlier
1060 data['outlier_rtc'] = dataOut.data_outlier
1053 1061
1054 1062 meta = {}
1055 1063
@@ -1059,7 +1067,7 class OutliersRTIPlot(Plot):
1059 1067 # self.data.normalize_heights()
1060 1068 self.x = self.data.times
1061 1069 self.y = self.data.yrange
1062 self.z = self.data['outlier']
1070 self.z = self.data['outlier_rtc']
1063 1071
1064 1072 #self.z = numpy.ma.masked_invalid(self.z)
1065 1073
@@ -1085,7 +1093,7 class OutliersRTIPlot(Plot):
1085 1093 cmap=self.cmaps[n]
1086 1094 )
1087 1095 if self.showprofile:
1088 ax.plot_profile = self.pf_axes[n].plot(data['outlier'][n], self.y)[0]
1096 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1089 1097 self.pf_axes[n].set_xlabel('')
1090 1098 else:
1091 1099 if self.zlimits is not None:
@@ -1097,84 +1105,83 class OutliersRTIPlot(Plot):
1097 1105 cmap=self.cmaps[n]
1098 1106 )
1099 1107 if self.showprofile:
1100 ax.plot_profile.set_data(data['outlier'][n], self.y)
1108 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1101 1109 self.pf_axes[n].set_xlabel('')
1102 1110
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)
1111 class NIncohIntRTIPlot(Plot):
1112 '''
1113 Plot for data_xxxx object
1114 '''
1115
1116 CODE = 'integrations_rtc' # Range Time Counts
1117 colormap = 'BuGn'
1118 plot_type = 'pcolorbuffer'
1119
1120 def setup(self):
1121 self.xaxis = 'time'
1122 self.ncols = 1
1123 self.nrows = self.data.shape('integrations_rtc')[0]
1124 self.nplots = self.nrows
1125 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1126
1127
1128 if not self.xlabel:
1129 self.xlabel = 'Time'
1130
1131 self.ylabel = 'Height [km]'
1132 if not self.titles:
1133 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1134
1135 def update(self, dataOut):
1136
1137 data = {}
1138 data['integrations_rtc'] = dataOut.nIncohInt
1139
1140 meta = {}
1141
1142 return data, meta
1143
1144 def plot(self):
1145 # self.data.normalize_heights()
1146 self.x = self.data.times
1147 self.y = self.data.yrange
1148 self.z = self.data['integrations_rtc']
1149
1150 #self.z = numpy.ma.masked_invalid(self.z)
1151
1152 if self.decimation is None:
1153 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1154 else:
1155 x, y, z = self.fill_gaps(*self.decimate())
1156
1157 for n, ax in enumerate(self.axes):
1158
1159 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1160 self.z[n])
1161 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1162 self.z[n])
1163 data = self.data[-1]
1164 if ax.firsttime:
1165 if self.zlimits is not None:
1166 self.zmin, self.zmax = self.zlimits[n]
1167
1168 ax.plt = ax.pcolormesh(x, y, z[n].T,
1169 vmin=self.zmin,
1170 vmax=self.zmax,
1171 cmap=self.cmaps[n]
1172 )
1173 if self.showprofile:
1174 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1175 self.pf_axes[n].set_xlabel('')
1176 else:
1177 if self.zlimits is not None:
1178 self.zmin, self.zmax = self.zlimits[n]
1179 ax.collections.remove(ax.collections[0])
1180 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1181 vmin=self.zmin,
1182 vmax=self.zmax,
1183 cmap=self.cmaps[n]
1184 )
1185 if self.showprofile:
1186 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1187 self.pf_axes[n].set_xlabel('')
Show file before | Show file after | Hide comments
@@ -465,7 +465,8 class HDFWriter(Operation):
465 465
466 466 def run(self, dataOut,**kwargs):
467 467
468 self.dataOut = dataOut.copy()
468 self.dataOut = dataOut
469
469 470 if not(self.isConfig):
470 471 self.setup(**kwargs)
471 472
Show file before | Show file after | Hide comments
@@ -59,7 +59,7 class ProcessingUnit(object):
59 59
60 60
61 61 if mybool:
62 #print("run jeje")
62 #print("run yeah")
63 63 self.run(**kwargs)
64 64 else:
65 65 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
@@ -210,7 +210,7 def MPDecorator(BaseClass):
210 210
211 211 def run(self):
212 212
213 while True:
213 while 1:
214 214
215 215 dataOut = self.queue.get()
216 216
Show file before | Show file after | Hide comments
@@ -101,7 +101,7 class ParametersProc(ProcessingUnit):
101 101 def run(self):
102 102 # print("run parameter proc")
103 103 #---------------------- Voltage Data ---------------------------
104
104 #print(self.dataIn.flagNoData)
105 105 if self.dataIn.type == "Voltage":
106 106
107 107 self.__updateObjFromInput()
@@ -132,10 +132,12 class ParametersProc(ProcessingUnit):
132 132 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
133 133 self.dataOut.data_spc = self.dataIn.data_spc
134 134 self.dataOut.data_cspc = self.dataIn.data_cspc
135 self.dataOut.data_outlier = self.dataIn.data_outlier
135 136 self.dataOut.nProfiles = self.dataIn.nProfiles
136 137 self.dataOut.nIncohInt = self.dataIn.nIncohInt
137 138 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
138 139 self.dataOut.ippFactor = self.dataIn.ippFactor
140 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
139 141 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
140 142 self.dataOut.ipp = self.dataIn.ipp
141 143 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
@@ -1472,9 +1474,9 class SpectralFitting(Operation):
1472 1474
1473 1475 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1474 1476
1475 print("run SpectralFitting")
1476 self.dataOut = dataOut.copy()
1477 print(self.dataOut.nIncohInt)
1477 #print("run SpectralFitting")
1478 self.dataOut = dataOut#.copy()
1479 #print(self.dataOut.nIncohInt)
1478 1480 if path != None:
1479 1481 sys.path.append(path)
1480 1482 self.dataOut.library = importlib.import_module(file)
Show file before | Show file after | Hide comments
@@ -21,7 +21,7 from schainpy.model.data import _noise
21 21
22 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 26 class SpectraProc(ProcessingUnit):
27 27
@@ -126,6 +126,11 class SpectraProc(ProcessingUnit):
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 try:
130 type = self.dataIn.type.decode("utf-8")
131 self.dataIn.type = type
132 except:
133 pass
129 134 if self.dataIn.type == "Spectra":
130 135
131 136 try:
@@ -209,9 +214,42 class SpectraProc(ProcessingUnit):
209 214 self.firstdatatime = None
210 215 self.profIndex = 0
211 216
217 elif self.dataIn.type == "Parameters":
218
219 self.dataOut.data_spc = self.dataIn.data_spc
220 self.dataOut.data_cspc = self.dataIn.data_cspc
221 self.dataOut.data_outlier = self.dataIn.data_outlier
222 self.dataOut.nProfiles = self.dataIn.nProfiles
223 self.dataOut.nIncohInt = self.dataIn.nIncohInt
224 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
225 self.dataOut.ippFactor = self.dataIn.ippFactor
226 self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
227 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
228 self.dataOut.ipp = self.dataIn.ipp
229 #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
230 #self.dataOut.spc_noise = self.dataIn.getNoise()
231 #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
232 # self.dataOut.normFactor = self.dataIn.normFactor
233 if hasattr(self.dataIn, 'channelList'):
234 self.dataOut.channelList = self.dataIn.channelList
235 if hasattr(self.dataIn, 'pairsList'):
236 self.dataOut.pairsList = self.dataIn.pairsList
237 self.dataOut.groupList = self.dataIn.pairsList
238
239 self.dataOut.flagNoData = False
240
241 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
242 self.dataOut.ChanDist = self.dataIn.ChanDist
243 else: self.dataOut.ChanDist = None
244
245 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
246 # self.dataOut.VelRange = self.dataIn.VelRange
247 #else: self.dataOut.VelRange = None
248
249
212 250
213 251 else:
214 raise ValueError("The type of input object '%s' is not valid".format(
252 raise ValueError("The type of input object {} is not valid".format(
215 253 self.dataIn.type))
216 254
217 255
@@ -497,14 +535,16 class removeDC(Operation):
497 535
498 536 return self.dataOut
499 537
500 class getNoise(Operation):
501 warnings = False
538 class getNoiseB(Operation):
539
540 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
502 541 def __init__(self):
503 542
504 543 Operation.__init__(self)
544 self.isConfig = False
545
546 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
505 547
506 def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
507 self.dataOut = dataOut
508 548 self.warnings = warnings
509 549 if minHei == None:
510 550 minHei = self.dataOut.heightList[0]
@@ -620,24 +660,111 class getNoise(Operation):
620 660 except:
621 661 maxIndexFFT = len( self.dataOut.getFreqRange(1))
622 662
663 self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
664 self.isConfig = True
665 if offset!=None:
666 self.offset = 10**(offset/10)
667 #print("config getNoise Done")
668
669 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
670 self.dataOut = dataOut
671
672 if not self.isConfig:
673 self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
623 674
624 675 self.dataOut.noise_estimation = None
625 676 noise = None
626 677 if self.dataOut.type == 'Voltage':
627 noise = self.dataOut.getNoise(ymin_index=minIndex, ymax_index=maxIndex)
678 noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
628 679 #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
629 680 elif self.dataOut.type == 'Spectra':
630 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
681
682 noise = numpy.zeros( self.dataOut.nChannels)
683 for channel in range( self.dataOut.nChannels):
684 norm = self.dataOut.max_nIncohInt/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
685 #print("norm nIncoh: ", norm )
686 daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
687 daux = numpy.multiply(daux, norm)
688 #print("offset: ", self.offset, 10*numpy.log10(self.offset))
689 #noise[channel] = self.getNoiseByMean(daux)/self.offset
690 noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
691
692 #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
631 693 else:
632 noise = self.dataOut.getNoise(xmin_index=minIndexFFT, xmax_index=maxIndexFFT, ymin_index=minIndex, ymax_index=maxIndex)
694 noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
633 695 self.dataOut.noise_estimation = noise.copy() # dataOut.noise
634 696 #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
635 697
698 #print(self.dataOut.flagNoData)
636 699 return self.dataOut
637 700
701 def getNoiseByMean(self,data):
702 #data debe estar ordenado
703 data = numpy.mean(data,axis=1)
704 sortdata = numpy.sort(data, axis=None)
705 #sortID=data.argsort()
706 #print(data.shape)
707
708 pnoise = None
709 j = 0
710
711 mean = numpy.mean(sortdata)
712 min = numpy.min(sortdata)
713 delta = mean - min
714 indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
715 #print(len(indexes))
716 if len(indexes)==0:
717 pnoise = numpy.mean(sortdata)
718 else:
719 j = indexes[0]
720 pnoise = numpy.mean(sortdata[0:j])
721
722 # from matplotlib import pyplot as plt
723 # plt.plot(sortdata)
724 # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r')
725 # plt.show()
726 #print("noise: ", 10*numpy.log10(pnoise))
727 return pnoise
728
729 def getNoiseByHS(self,data, navg):
730 #data debe estar ordenado
731 #data = numpy.mean(data,axis=1)
732 sortdata = numpy.sort(data, axis=None)
733
734 lenOfData = len(sortdata)
735 nums_min = lenOfData*0.05
736
737 if nums_min <= 5:
738
739 nums_min = 5
740
741 sump = 0.
742 sumq = 0.
743
744 j = 0
745 cont = 1
746
747 while((cont == 1)and(j < lenOfData)):
748
749 sump += sortdata[j]
750 sumq += sortdata[j]**2
751 #sumq -= sump**2
752 if j > nums_min:
753 rtest = float(j)/(j-1) + 1.0/0.1
754 #if ((sumq*j) > (sump**2)):
755 if ((sumq*j) > (rtest*sump**2)):
756 j = j - 1
757 sump = sump - sortdata[j]
758 sumq = sumq - sortdata[j]**2
759 cont = 0
760
761 j += 1
762
763 lnoise = sump / j
764
765 return lnoise
638 766
639 767
640 # import matplotlib.pyplot as plt
641 768
642 769 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
643 770 z = (x - a1) / a2
@@ -779,7 +906,7 class CleanRayleigh(Operation):
779 906
780 907 #index = tini.tm_hour*12+tini.tm_min/5
781 908 '''
782 REVISAR
909 #REVISAR
783 910 '''
784 911 # jspc = jspc/self.nFFTPoints/self.normFactor
785 912 # jcspc = jcspc/self.nFFTPoints/self.normFactor
@@ -1396,7 +1523,7 class IntegrationFaradaySpectra(Operation):
1396 1523 data_dc = self.__buffer_dc
1397 1524 #(CH, HEIGH)
1398 1525 self.maxProfilesInt = self.__profIndex
1399 n = self.__profIndex - self.dataOutliers
1526 n = self.__profIndex - self.dataOutliers # n becomes a matrix
1400 1527
1401 1528 self.__buffer_spc = []
1402 1529 self.__buffer_cspc = []
@@ -1455,11 +1582,11 class IntegrationFaradaySpectra(Operation):
1455 1582 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1456 1583
1457 1584 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1458 self.dataOut = dataOut.copy()
1585 self.dataOut = dataOut
1459 1586 if n == 1:
1460 1587 return self.dataOut
1461 1588
1462 self.dataOut.flagNoData = True
1589
1463 1590 if self.dataOut.nChannels == 1:
1464 1591 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1465 1592 #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc)
@@ -1483,7 +1610,7 class IntegrationFaradaySpectra(Operation):
1483 1610 self.dataOut.dataLag_spc,
1484 1611 self.dataOut.dataLag_cspc,
1485 1612 self.dataOut.dataLag_dc)
1486
1613 self.dataOut.flagNoData = True
1487 1614 if self.__dataReady:
1488 1615
1489 1616 if not self.ByLags:
@@ -1513,6 +1640,7 class IntegrationFaradaySpectra(Operation):
1513 1640 self.dataOut.utctime = avgdatatime
1514 1641 self.dataOut.flagNoData = False
1515 1642 #print("Faraday Integration DONE...")
1643 #print(self.dataOut.flagNoData)
1516 1644 return self.dataOut
1517 1645
1518 1646 class removeInterference(Operation):
@@ -1910,7 +2038,7 class IncohInt(Operation):
1910 2038 self.incohInt += dataOut.nIncohInt
1911 2039 self.nOutliers += dataOut.data_outlier
1912 2040 if self.__dataReady:
1913
2041 #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
1914 2042 dataOut.data_spc = avgdata_spc
1915 2043 dataOut.data_cspc = avgdata_cspc
1916 2044 dataOut.data_dc = avgdata_dc
@@ -1918,7 +2046,7 class IncohInt(Operation):
1918 2046 dataOut.data_outlier = self.nOutliers
1919 2047 dataOut.utctime = avgdatatime
1920 2048 dataOut.flagNoData = False
1921 dataOut.max_nIncohInt = self.__profIndex
2049 dataOut.max_nIncohInt += self.__profIndex
1922 2050 self.incohInt = 0
1923 2051 self.nOutliers = 0
1924 2052 self.__profIndex = 0
Show file before | Show file after | Hide comments
@@ -6,7 +6,7 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 7 from schainpy.model.io.utils import getHei_index
8 8 from time import time
9 import datetime
9 #import datetime
10 10 import numpy
11 11 #import copy
12 12 from schainpy.model.data import _noise
@@ -1742,6 +1742,7 class SSheightProfiles(Operation):
1742 1742 dataOut.step = self.step
1743 1743 #print(numpy.shape(dataOut.data))
1744 1744 #exit(1)
1745 #print("new data shape and time:", dataOut.data.shape, dataOut.utctime)
1745 1746
1746 1747 return dataOut
1747 1748 ################################################################################3############################3
@@ -1754,14 +1755,17 class SSheightProfiles2(Operation):
1754 1755 Procesa por perfiles y por bloques
1755 1756 '''
1756 1757
1757 step = None
1758 nsamples = None
1758
1759 1759 bufferShape = None
1760 1760 profileShape = None
1761 1761 sshProfiles = None
1762 1762 profileIndex = None
1763 deltaHeight = None
1764 init_range = None
1763 #nsamples = None
1764 #step = None
1765 #deltaHeight = None
1766 #init_range = None
1767 __slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels',
1768 '__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights')
1765 1769
1766 1770 def __init__(self, **kwargs):
1767 1771
@@ -1786,12 +1790,12 class SSheightProfiles2(Operation):
1786 1790 self.deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1787 1791 self.init_range = dataOut.heightList[0]
1788 1792 #numberProfile = self.nsamples
1789 self.numberSamples = (self.__nHeis - self.nsamples)/self.step
1793 numberSamples = (self.__nHeis - self.nsamples)/self.step
1790 1794
1791 self.new_nHeights = self.numberSamples
1795 self.new_nHeights = numberSamples
1792 1796
1793 self.bufferShape = int(self.__nChannels), int(self.numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
1794 self.profileShape = int(self.__nChannels), int(self.nsamples), int(self.numberSamples) # nchannels, nprofiles, nsamples
1797 self.bufferShape = int(self.__nChannels), int(numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
1798 self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples) # nchannels, nprofiles, nsamples
1795 1799
1796 1800 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1797 1801 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
@@ -1878,7 +1882,8 class SSheightProfiles2(Operation):
1878 1882 dataOut.flagDataAsBlock = True
1879 1883
1880 1884 dataBlock = None
1881 #print("new data shape:", dataOut.data.shape)
1885
1886 #print("new data shape:", dataOut.data.shape, dataOut.utctime)
1882 1887
1883 1888 return dataOut
1884 1889
@@ -1893,27 +1898,28 class removeProfileByFaradayHS(Operation):
1893 1898 '''
1894 1899
1895 1900 '''
1896 isConfig = False
1897 n = None
1898
1899 __profIndex = 0
1901 #isConfig = False
1902 #n = None
1900 1903
1901 __dataReady = False
1904 #__dataReady = False
1902 1905 __buffer_data = []
1903 1906 __buffer_times = []
1904 __initime = None
1905 __count_exec = 0
1906 __profIndex = 0
1907 #__initime = None
1908 #__count_exec = 0
1909 #__profIndex = 0
1907 1910 buffer = None
1908 lenProfileOut = 1
1911 #lenProfileOut = 1
1912
1913 #init_prof = 0
1914 #end_prof = 0
1909 1915
1910 init_prof = 0
1911 end_prof = 0
1912 n_profiles = 0
1913 first_utcBlock = None
1916 #first_utcBlock = None
1914 1917 outliers_IDs_list = []
1915 __dh = 0
1918 #__dh = 0
1916 1919
1920 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
1921 '__dh','first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
1922 '__count_exec','__initime','__dataReady','__ipp')
1917 1923 def __init__(self, **kwargs):
1918 1924
1919 1925 Operation.__init__(self, **kwargs)
@@ -1932,12 +1938,13 class removeProfileByFaradayHS(Operation):
1932 1938 self.thHistOutlier = thHistOutlier
1933 1939 self.__profIndex = 0
1934 1940 self.buffer = None
1935 self.lenProfileOut = 1
1941 self._ipp = dataOut.ippSeconds
1936 1942 self.n_prof_released = 0
1937 1943 self.heightList = dataOut.heightList
1938 1944 self.init_prof = 0
1939 1945 self.end_prof = 0
1940 self.n_profiles = 0
1946 self.__count_exec = 0
1947 self.__profIndex = 0
1941 1948 self.first_utcBlock = None
1942 1949 self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
1943 1950 minHei = minHei
@@ -1948,6 +1955,8 class removeProfileByFaradayHS(Operation):
1948 1955 maxHei = dataOut.heightList[-1]
1949 1956 self.minHei_idx,self.maxHei_idx = getHei_index(minHei, maxHei, dataOut.heightList)
1950 1957
1958 self.nChannels = dataOut.nChannels
1959 self.nHeights = dataOut.nHeights
1951 1960
1952 1961 def filterSatsProfiles(self):
1953 1962 data = self.__buffer_data
@@ -2124,7 +2133,7 class removeProfileByFaradayHS(Operation):
2124 2133 def releaseBlock(self):
2125 2134
2126 2135 if self.n % self.lenProfileOut != 0:
2127 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2136 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2128 2137 return None
2129 2138
2130 2139 data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
@@ -2139,9 +2148,9 class removeProfileByFaradayHS(Operation):
2139 2148 return data
2140 2149
2141 2150 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None):
2142 #print("run op buffer 2D")
2143 self.nChannels = dataOut.nChannels
2144 self.nHeights = dataOut.nHeights
2151 #print("run op buffer 2D",dataOut.ippSeconds)
2152 # self.nChannels = dataOut.nChannels
2153 # self.nHeights = dataOut.nHeights
2145 2154
2146 2155 if not self.isConfig:
2147 2156 #print("init p idx: ", dataOut.profileIndex )
@@ -2158,9 +2167,7 class removeProfileByFaradayHS(Operation):
2158 2167 self.lenProfileOut = nProfilesOut
2159 2168 dataOut.flagNoData = False
2160 2169 #print("tp 2 ",dataOut.data.shape)
2161 #(ch, self.n_profiles, nh) = self.buffer.shape
2162 #print("tp 3 ",self.dataOut.data.shape)
2163 #print("rel: ",self.buffer[:,-1,:])
2170
2164 2171 self.init_prof = 0
2165 2172 self.end_prof = self.lenProfileOut
2166 2173
@@ -2174,12 +2181,13 class removeProfileByFaradayHS(Operation):
2174 2181 if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
2175 2182 #print("omitting: ", self.n_prof_released)
2176 2183 dataOut.flagNoData = True
2177
2178 dataOut.utctime = self.first_utcBlock + self.init_prof*dataOut.ippSeconds
2184 dataOut.ippSeconds = self._ipp
2185 dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
2186 # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
2179 2187 #dataOut.data = self.releaseBlock()
2180 2188 #########################################################3
2181 2189 if self.n % self.lenProfileOut != 0:
2182 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2190 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
2183 2191 return None
2184 2192
2185 2193 dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:] #ch, prof, alt
General Comments 0
You need to be logged in to leave comments. Login now