schain Commit - r1475:0c13fcec0e07 · Jicamarca Repository

Hildebrand_S en C, clear sats

joabAM -

r1475:0c13fcec0e07

parent child

Context file:

r1475:0c13fcec0e07

Collapse all files

schainc/_noise.c +11 -11

               if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
                   return NULL;
               }
               data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
               if (data_array == NULL) {
                   Py_XDECREF(data_array);
                   Py_XDECREF(data_obj);
               // return PyLong_FromLong(lnoise);
               return PyFloat_FromDouble(lnoise);
             }
-            /*
             static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
               double navg;
-              double th;
               PyObject *data_obj, *data_array;
-              if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg, &th)) {
+              if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
                   return NULL;
               }
               data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
               if (data_array == NULL) {
                   Py_XDECREF(data_array);
                   Py_XDECREF(data_obj);
               }
               double *sortdata = (double*)PyArray_DATA(data_array);
               int lenOfData = (int)PyArray_SIZE(data_array) ;
-              double nums_min = lenOfData*th;
+              double nums_min = lenOfData*0.75;
               if (nums_min <= 5) nums_min = 5;
               double sump = 0;
               double sumq = 0;
               Py_DECREF(data_array);
-              // return PyLong_FromLong(lnoise);
+              return PyLong_FromLong(j);
-              return PyFloat_FromDouble(j,sortID);
             }
-            */
             static PyMethodDef noiseMethods[] = {
               { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
+              { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Get index for satellite cleaning" },
               { NULL, NULL, 0, NULL }
             };

schainpy/model/data/jrodata.py +4 -4

                 data = None
                 nmodes = None
                 metadata_list = ['heightList', 'timeZone', 'type']
-                codeList = None
+                codeList = []
-                azimuthList = None
+                azimuthList = []
-                elevationList = None
+                elevationList = []
                 def __str__(self):
                     return noise
                 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     if self.noise_estimation is not None:
                         # this was estimated by getNoise Operation defined in jroproc_spectra.py
                         return self.noise_estimation

schainpy/model/graphics/jroplot_parameters.py +4 -1

                     data = {
                         'pow': 10*numpy.log10(dataOut.data_pow)
                     }
-                    data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
+                    try:
+                        data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
+                    except:
+                        pass
                     return data, {}
             class SpectralWidthPlot(RTIPlot):

schainpy/model/graphics/jroplot_spectra.py +50 -12

                 plot_type = 'pcolor'
                 buffering = False
                 channelList = []
+                elevationList = []
+                azimuthList = []
                 def setup(self):
                 def update_list(self,dataOut):
                     if len(self.channelList) == 0:
                         self.channelList = dataOut.channelList
+                    if len(self.elevationList) == 0:
+                        self.elevationList = dataOut.elevationList
+                    if len(self.azimuthList) == 0:
+                        self.azimuthList = dataOut.azimuthList
                 def update(self, dataOut):
                                 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
                             if self.CODE == 'spc_moments':
                                 ax.plt_mean.set_data(mean, y)
-                        self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
+                        if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
+                            self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
+                        else:
+                            self.titles.append('CH {}:  {:3.2f}dB'.format(self.channelList[n], noise))
             class CrossSpectraPlot(Plot):
                 plot_type = 'pcolorbuffer'
                 titles = None
                 channelList = []
+                elevationList = []
+                azimuthList = []
                 def setup(self):
                     self.xaxis = 'time'
                 def update_list(self,dataOut):
-                    self.channelList = dataOut.channelList
+                    if len(self.channelList) == 0:
+                        self.channelList = dataOut.channelList
+                    if len(self.elevationList) == 0:
+                        self.elevationList = dataOut.elevationList
+                    if len(self.azimuthList) == 0:
+                        self.azimuthList = dataOut.azimuthList
                 def update(self, dataOut):
                     try:
                         if self.channelList != None:
-                            self.titles = ['{} Channel {}'.format(
+                            if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
-                                self.CODE.upper(), x) for x in self.channelList]
+                                self.titles = ['{} Channel {} ({:2.1f} Elev,  {:2.1f} Azth)'.format(
+                                    self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
+                            else:
+                                self.titles = ['{} Channel {}'.format(
+                                    self.CODE.upper(), x) for x in self.channelList]
                     except:
                         if self.channelList.any() != None:
                             self.titles = ['{} Channel {}'.format(
                                 self.CODE.upper(), x) for x in self.channelList]
                     if self.decimation is None:
                         x, y, z = self.fill_gaps(self.x, self.y, self.z)
                     else:
                             labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
                             self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
                         else:
                             for i, line in enumerate(ax.plt):
                                 line.set_data(x, z[n, :, index[i]])
                         self.titles.append('CH {}'.format(self.channelList[n]))
                     plt.suptitle(self.maintitle,  fontsize=10)
                 plot_type = 'pcolorbuffer'
                 titles = None
                 channelList = []
+                elevationList = []
+                azimuthList = []
                 def setup(self):
                     self.xaxis = 'time'
                         self.CODE.upper(), x) for x in range(self.nplots)]
                 def update_list(self,dataOut):
+                    if len(self.channelList) == 0:
-                    self.channelList = dataOut.channelList
+                        self.channelList = dataOut.channelList
+                    if len(self.elevationList) == 0:
+                        self.elevationList = dataOut.elevationList
+                    if len(self.azimuthList) == 0:
+                        self.azimuthList = dataOut.azimuthList
                 def update(self, dataOut):
                     if len(self.channelList) == 0:
                     n0 = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     (nch, nff, nh) = dataOut.data_spc.shape
+                    #print(nch, nff, nh)
+                    if nch != 1:
+                        aux = []
+                        for c in self.channelList:
+                            aux.append(n0[c])
+                        n0 = numpy.asarray(aux)
                     noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
+                    #print(dataOut.elevationList, dataOut.azimuthList)
+                    #print(dataOut.channelList)
                     data['noiseless_rti'] = dataOut.getPower() - noise
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     return data, meta
                     try:
                         if self.channelList != None:
-                            self.titles = ['{} Channel {}'.format(
+                            if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
-                                self.CODE.upper(), x) for x in self.channelList]
+                                self.titles = ['{} Channel {} ({:2.1f} Elev,  {:2.1f} Azth)'.format(
+                                    self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
+                            else:
+                                self.titles = ['{} Channel {}'.format(
+                                    self.CODE.upper(), x) for x in self.channelList]
                     except:
                         if self.channelList.any() != None:
                             self.titles = ['{} Channel {}'.format(
                                 self.CODE.upper(), x) for x in self.channelList]
                     if self.decimation is None:

schainpy/model/proc/jroproc_spectra.py +16 -7

             from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
             from schainpy.model.data.jrodata import Spectra
             from schainpy.model.data.jrodata import hildebrand_sekhon
+            from schainpy.model.data import _noise
             from schainpy.utils import log
             from scipy.optimize import curve_fit
                 n = None
                 minHei_ind = None
                 maxHei_ind = None
+                avg = 1.0
                 factor = 0.0
                 def __init__(self):
                     return
-                def hildebrand_sekhon_Integration(self,data,navg, factor):
+                def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
+                    #data debe estar ordenado
-                    sortdata = numpy.sort(data, axis=None)
+                    #sortdata = numpy.sort(data, axis=None)
-                    sortID=data.argsort()
+                    #sortID=data.argsort()
                     lenOfData = len(sortdata)
                     nums_min = lenOfData*factor
                     if nums_min <= 5:
                         j += 1
                     #lnoise = sump / j
                     #print("H S done")
-                    return j,sortID
+                    #return j,sortID
+                    return j
                 def pushData(self):
                     """
                                 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
                                 #     continue
                                 buffer=buffer1[:,j]
-                                index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
+                                sortdata = numpy.sort(buffer, axis=None)
+                                sortID=buffer.argsort()
+                                index = _noise.hildebrand_sekhon2(sortdata,self.navg)
+                                #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
                                 indexes.append(index)
                                 #sortIDs.append(sortID)
                             if indexmin != buffer1.shape[0]:
                                 if self.nChannels > 1:
                                     cspc_outliers_exist= True
-                                    print("outliers cpsc")
+                                    #print("outliers cspc")
                                 ###sortdata=numpy.sort(buffer1,axis=0)
                                 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
                                 lt=outliers_IDs

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages