##// END OF EJS Templates
schain
RSS
Hildebrand_S en C, clear sats
joabAM -
r1475:0c13fcec0e07
parent child
Show More
Side by Side Unified
Context file:
r1475:0c13fcec0e07 Loading diff...:
Show file before | Show file after | Hide comments
@@ -10,9 +10,9 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
10 10 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
11 11 return NULL;
12 12 }
13
13
14 14 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
15
15
16 16 if (data_array == NULL) {
17 17 Py_XDECREF(data_array);
18 18 Py_XDECREF(data_obj);
@@ -49,18 +49,17 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
49 49 // return PyLong_FromLong(lnoise);
50 50 return PyFloat_FromDouble(lnoise);
51 51 }
52 /*
52
53 53 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
54 54 double navg;
55 double th;
56 55 PyObject *data_obj, *data_array;
57 56
58 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg, &th)) {
57 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
59 58 return NULL;
60 59 }
61
60
62 61 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
63
62
64 63 if (data_array == NULL) {
65 64 Py_XDECREF(data_array);
66 65 Py_XDECREF(data_obj);
@@ -68,7 +67,7 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
68 67 }
69 68 double *sortdata = (double*)PyArray_DATA(data_array);
70 69 int lenOfData = (int)PyArray_SIZE(data_array) ;
71 double nums_min = lenOfData*th;
70 double nums_min = lenOfData*0.75;
72 71 if (nums_min <= 5) nums_min = 5;
73 72 double sump = 0;
74 73 double sumq = 0;
@@ -94,13 +93,14 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
94 93
95 94 Py_DECREF(data_array);
96 95
97 // return PyLong_FromLong(lnoise);
98 return PyFloat_FromDouble(j,sortID);
96 return PyLong_FromLong(j);
97
99 98 }
100 */
99
101 100
102 101 static PyMethodDef noiseMethods[] = {
103 102 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
103 { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Get index for satellite cleaning" },
104 104 { NULL, NULL, 0, NULL }
105 105 };
106 106
Show file before | Show file after | Hide comments
@@ -192,9 +192,9 class JROData(GenericData):
192 192 data = None
193 193 nmodes = None
194 194 metadata_list = ['heightList', 'timeZone', 'type']
195 codeList = None
196 azimuthList = None
197 elevationList = None
195 codeList = []
196 azimuthList = []
197 elevationList = []
198 198
199 199 def __str__(self):
200 200
@@ -489,7 +489,7 class Spectra(JROData):
489 489 return noise
490 490
491 491 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
492
492
493 493 if self.noise_estimation is not None:
494 494 # this was estimated by getNoise Operation defined in jroproc_spectra.py
495 495 return self.noise_estimation
Show file before | Show file after | Hide comments
@@ -87,7 +87,10 class PowerPlot(RTIPlot):
87 87 data = {
88 88 'pow': 10*numpy.log10(dataOut.data_pow)
89 89 }
90 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
90 try:
91 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
92 except:
93 pass
91 94 return data, {}
92 95
93 96 class SpectralWidthPlot(RTIPlot):
Show file before | Show file after | Hide comments
@@ -23,6 +23,8 class SpectraPlot(Plot):
23 23 plot_type = 'pcolor'
24 24 buffering = False
25 25 channelList = []
26 elevationList = []
27 azimuthList = []
26 28
27 29 def setup(self):
28 30
@@ -43,6 +45,10 class SpectraPlot(Plot):
43 45 def update_list(self,dataOut):
44 46 if len(self.channelList) == 0:
45 47 self.channelList = dataOut.channelList
48 if len(self.elevationList) == 0:
49 self.elevationList = dataOut.elevationList
50 if len(self.azimuthList) == 0:
51 self.azimuthList = dataOut.azimuthList
46 52
47 53 def update(self, dataOut):
48 54
@@ -110,7 +116,10 class SpectraPlot(Plot):
110 116 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
111 117 if self.CODE == 'spc_moments':
112 118 ax.plt_mean.set_data(mean, y)
113 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
119 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
120 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
121 else:
122 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
114 123
115 124
116 125 class CrossSpectraPlot(Plot):
@@ -226,6 +235,8 class RTIPlot(Plot):
226 235 plot_type = 'pcolorbuffer'
227 236 titles = None
228 237 channelList = []
238 elevationList = []
239 azimuthList = []
229 240
230 241 def setup(self):
231 242 self.xaxis = 'time'
@@ -242,7 +253,12 class RTIPlot(Plot):
242 253
243 254 def update_list(self,dataOut):
244 255
245 self.channelList = dataOut.channelList
256 if len(self.channelList) == 0:
257 self.channelList = dataOut.channelList
258 if len(self.elevationList) == 0:
259 self.elevationList = dataOut.elevationList
260 if len(self.azimuthList) == 0:
261 self.azimuthList = dataOut.azimuthList
246 262
247 263
248 264 def update(self, dataOut):
@@ -265,12 +281,18 class RTIPlot(Plot):
265 281
266 282 try:
267 283 if self.channelList != None:
268 self.titles = ['{} Channel {}'.format(
269 self.CODE.upper(), x) for x in self.channelList]
284 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
285 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
286 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
287 else:
288 self.titles = ['{} Channel {}'.format(
289 self.CODE.upper(), x) for x in self.channelList]
270 290 except:
271 291 if self.channelList.any() != None:
292
272 293 self.titles = ['{} Channel {}'.format(
273 294 self.CODE.upper(), x) for x in self.channelList]
295
274 296 if self.decimation is None:
275 297 x, y, z = self.fill_gaps(self.x, self.y, self.z)
276 298 else:
@@ -532,7 +554,7 class SpectraCutPlot(Plot):
532 554 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
533 555 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
534 556 else:
535 for i, line in enumerate(ax.plt):
557 for i, line in enumerate(ax.plt):
536 558 line.set_data(x, z[n, :, index[i]])
537 559 self.titles.append('CH {}'.format(self.channelList[n]))
538 560 plt.suptitle(self.maintitle, fontsize=10)
@@ -879,6 +901,8 class NoiselessRTIPlot(Plot):
879 901 plot_type = 'pcolorbuffer'
880 902 titles = None
881 903 channelList = []
904 elevationList = []
905 azimuthList = []
882 906
883 907 def setup(self):
884 908 self.xaxis = 'time'
@@ -894,9 +918,12 class NoiselessRTIPlot(Plot):
894 918 self.CODE.upper(), x) for x in range(self.nplots)]
895 919
896 920 def update_list(self,dataOut):
897
898 self.channelList = dataOut.channelList
899
921 if len(self.channelList) == 0:
922 self.channelList = dataOut.channelList
923 if len(self.elevationList) == 0:
924 self.elevationList = dataOut.elevationList
925 if len(self.azimuthList) == 0:
926 self.azimuthList = dataOut.azimuthList
900 927
901 928 def update(self, dataOut):
902 929 if len(self.channelList) == 0:
@@ -906,9 +933,15 class NoiselessRTIPlot(Plot):
906 933
907 934 n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
908 935 (nch, nff, nh) = dataOut.data_spc.shape
936 #print(nch, nff, nh)
937 if nch != 1:
938 aux = []
939 for c in self.channelList:
940 aux.append(n0[c])
941 n0 = numpy.asarray(aux)
909 942 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
910
911
943 #print(dataOut.elevationList, dataOut.azimuthList)
944 #print(dataOut.channelList)
912 945 data['noiseless_rti'] = dataOut.getPower() - noise
913 946 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
914 947 return data, meta
@@ -923,10 +956,15 class NoiselessRTIPlot(Plot):
923 956
924 957 try:
925 958 if self.channelList != None:
926 self.titles = ['{} Channel {}'.format(
927 self.CODE.upper(), x) for x in self.channelList]
959 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
960 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
961 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
962 else:
963 self.titles = ['{} Channel {}'.format(
964 self.CODE.upper(), x) for x in self.channelList]
928 965 except:
929 966 if self.channelList.any() != None:
967
930 968 self.titles = ['{} Channel {}'.format(
931 969 self.CODE.upper(), x) for x in self.channelList]
932 970 if self.decimation is None:
Show file before | Show file after | Hide comments
@@ -17,6 +17,8 import math
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 from schainpy.model.data import _noise
21
20 22 from schainpy.utils import log
21 23
22 24 from scipy.optimize import curve_fit
@@ -1092,6 +1094,7 class IntegrationFaradaySpectra(Operation):
1092 1094 n = None
1093 1095 minHei_ind = None
1094 1096 maxHei_ind = None
1097 avg = 1.0
1095 1098 factor = 0.0
1096 1099
1097 1100 def __init__(self):
@@ -1184,10 +1187,10 class IntegrationFaradaySpectra(Operation):
1184 1187
1185 1188 return
1186 1189
1187 def hildebrand_sekhon_Integration(self,data,navg, factor):
1188
1189 sortdata = numpy.sort(data, axis=None)
1190 sortID=data.argsort()
1190 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
1191 #data debe estar ordenado
1192 #sortdata = numpy.sort(data, axis=None)
1193 #sortID=data.argsort()
1191 1194 lenOfData = len(sortdata)
1192 1195 nums_min = lenOfData*factor
1193 1196 if nums_min <= 5:
@@ -1209,7 +1212,9 class IntegrationFaradaySpectra(Operation):
1209 1212 j += 1
1210 1213 #lnoise = sump / j
1211 1214 #print("H S done")
1212 return j,sortID
1215 #return j,sortID
1216 return j
1217
1213 1218
1214 1219 def pushData(self):
1215 1220 """
@@ -1256,7 +1261,11 class IntegrationFaradaySpectra(Operation):
1256 1261 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1257 1262 # continue
1258 1263 buffer=buffer1[:,j]
1259 index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1264 sortdata = numpy.sort(buffer, axis=None)
1265 sortID=buffer.argsort()
1266 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
1267
1268 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
1260 1269
1261 1270 indexes.append(index)
1262 1271 #sortIDs.append(sortID)
@@ -1272,7 +1281,7 class IntegrationFaradaySpectra(Operation):
1272 1281 if indexmin != buffer1.shape[0]:
1273 1282 if self.nChannels > 1:
1274 1283 cspc_outliers_exist= True
1275 print("outliers cpsc")
1284 #print("outliers cspc")
1276 1285 ###sortdata=numpy.sort(buffer1,axis=0)
1277 1286 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1278 1287 lt=outliers_IDs
General Comments 0
You need to be logged in to leave comments. Login now