##// END OF EJS Templates
schain
RSS
omisión perfies para ISR, Bloques a perfiles, heildebrand en C, y más...
joabAM -
r1506:7d81a1c691b1
parent child
Show More
Side by Side Unified
Context file:
r1506:7d81a1c691b1 Loading diff...:
Show file before | Show file after | Hide comments
@@ -57,7 +57,6 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
57 57 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
58 58 return NULL;
59 59 }
60
61 60 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
62 61
63 62 if (data_array == NULL) {
@@ -89,18 +88,15 static PyObject *hildebrand_sekhon2(PyObject *self, PyObject *args) {
89 88 j = j + 1;
90 89 }
91 90
92 //double lnoise = sump / j;
93
94 91 Py_DECREF(data_array);
95 92
96 93 return PyLong_FromLong(j);
97
98 94 }
99 95
100 96
101 97 static PyMethodDef noiseMethods[] = {
102 98 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
103 { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Get index for satellite cleaning" },
99 { "hildebrand_sekhon2", hildebrand_sekhon2, METH_VARARGS, "Variation for satellite cleaning" },
104 100 { NULL, NULL, 0, NULL }
105 101 };
106 102
Show file before | Show file after | Hide comments
@@ -266,10 +266,8 class ProcUnitConf(ConfBase):
266 266 def run(self):
267 267 '''
268 268 '''
269
270 269 return self.object.call(**self.getKwargs())
271 270
272
273 271 class ReadUnitConf(ProcUnitConf):
274 272
275 273 ELEMENTNAME = 'ReadUnit'
@@ -633,18 +631,47 class Project(Process):
633 631
634 632 err = False
635 633 n = len(self.configurations)
636
634 flag_no_read = False
635 nProc_noRead = 0
637 636 while not err:
638 637 #print("STAR")
638 n_proc = 0
639 639 for conf in self.getUnits():
640 640 #print("CONF: ",conf)
641 ok = conf.run()
641 #print( n_proc, nProc_noRead,flag_no_read)
642 if flag_no_read:
643 if n_proc >= nProc_noRead:
644 ok = conf.run()
645 else:
646 #print("break")
647 #print(ok, n_proc, nProc_noRead, n, flag_no_read)
648 n_proc += 1
649 continue
650 else:
651 ok = conf.run()
652
653 n_proc += 1
654 #print(ok, n_proc, nProc_noRead, n, flag_no_read)
655
642 656 if ok == 'Error':
643 657 #self.removeProcUnit(conf.id) #remove proc Unit
644 658 n -= 1
645 659 continue
660
661 elif ok == 'no_Read' and (not flag_no_read):
662 nProc_noRead = n_proc - 1
663 flag_no_read = True
664 print("No read")
665 continue
666 elif ok == 'new_Read':
667 nProc_noRead = 0
668 flag_no_read = False
669 print("read again")
670 continue
646 671 elif not ok:
672 #print("not ok",ok)
647 673 break
674
648 675 if n == 0:
649 676 err = True
650 677
Show file before | Show file after | Hide comments
@@ -159,6 +159,9 class GenericData(object):
159 159
160 160 class JROData(GenericData):
161 161
162 useInputBuffer = False
163 buffer_empty = True
164
162 165 systemHeaderObj = SystemHeader()
163 166 radarControllerHeaderObj = RadarControllerHeader()
164 167 type = None
Show file before | Show file after | Hide comments
@@ -434,6 +434,7 class Plot(Operation):
434 434 ax.firsttime = False
435 435 if self.grid:
436 436 ax.grid(True)
437
437 438 if not self.polar:
438 439 ax.set_title('{} {} {}'.format(
439 440 self.titles[n],
@@ -442,6 +443,7 class Plot(Operation):
442 443 self.time_label),
443 444 size=8)
444 445 else:
446
445 447 ax.set_title('{}'.format(self.titles[n]), size=8)
446 448 ax.set_ylim(0, 90)
447 449 ax.set_yticks(numpy.arange(0, 90, 20))
@@ -480,6 +482,7 class Plot(Operation):
480 482
481 483 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
482 484 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
485
483 486 fig.canvas.draw()
484 487 if self.show:
485 488 fig.show()
Show file before | Show file after | Hide comments
@@ -11,7 +11,7 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13 from itertools import combinations
14
14 from matplotlib.ticker import LinearLocator
15 15
16 16 class SpectraPlot(Plot):
17 17 '''
@@ -86,7 +86,7 class SpectraPlot(Plot):
86 86
87 87 data = self.data[-1]
88 88 z = data['spc']
89
89 print(z.shape, x.shape, y.shape)
90 90 for n, ax in enumerate(self.axes):
91 91 noise = data['noise'][n]
92 92 if self.CODE == 'spc_moments':
@@ -494,11 +494,13 class SpectraCutPlot(Plot):
494 494 self.height = 4.8 * self.nrows
495 495 self.ylabel = 'Power [dB]'
496 496 self.colorbar = False
497 self.plots_adjust.update({'left':0.15, 'hspace':0.3, 'right': 0.85, 'bottom':0.08})
497 self.plots_adjust.update({'left':0.05, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
498 498
499 499 if len(self.selectedHeightsList) > 0:
500 500 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
501 501
502
503
502 504 def update(self, dataOut):
503 505 if len(self.channelList) == 0:
504 506 self.channelList = dataOut.channelList
@@ -544,21 +546,33 class SpectraCutPlot(Plot):
544 546 else:
545 547 index = numpy.arange(0, len(y), int((len(y))/9))
546 548 #print("inde x ", index, self.axes)
549
547 550 for n, ax in enumerate(self.axes):
551
548 552 if ax.firsttime:
553
554
549 555 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
550 556 self.xmin = self.xmin if self.xmin else -self.xmax
551 557 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
552 558 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
559
560
553 561 ax.plt = ax.plot(x, z[n, :, index].T)
554 562 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
555 563 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
564 ax.minorticks_on()
565 ax.grid(which='major', axis='both')
566 ax.grid(which='minor', axis='x')
556 567 else:
557 568 for i, line in enumerate(ax.plt):
558 569 line.set_data(x, z[n, :, index[i]])
570
571
559 572 self.titles.append('CH {}'.format(self.channelList[n]))
560 573 plt.suptitle(self.maintitle, fontsize=10)
561 574
575
562 576 class BeaconPhase(Plot):
563 577
564 578 __isConfig = None
@@ -937,7 +951,7 class NoiselessRTIPlot(Plot):
937 951 if nch != 1:
938 952 aux = []
939 953 for c in self.channelList:
940 aux.append(n0[c])
954 aux.append(n0[c])
941 955 n0 = numpy.asarray(aux)
942 956 noise = numpy.repeat(n0,nh, axis=0).reshape((nch,nh))
943 957 #print(dataOut.elevationList, dataOut.azimuthList)
Show file before | Show file after | Hide comments
@@ -4,6 +4,7 Utilities for IO modules
4 4
5 5 import os
6 6 from datetime import datetime
7 import numpy
7 8
8 9 def folder_in_range(folder, start_date, end_date, pattern):
9 10 """
@@ -22,3 +23,29 def folder_in_range(folder, start_date, end_date, pattern):
22 23 except:
23 24 raise ValueError('Folder {} does not match {} format'.format(folder, pattern))
24 25 return start_date <= dt.date() <= end_date
26
27
28 def getHei_index( minHei, maxHei, heightList):
29 if (minHei < heightList[0]):
30 minHei = heightList[0]
31
32 if (maxHei > heightList[-1]):
33 maxHei = heightList[-1]
34
35 minIndex = 0
36 maxIndex = 0
37 heights = numpy.asarray(heightList)
38
39 inda = numpy.where(heights >= minHei)
40 indb = numpy.where(heights <= maxHei)
41
42 try:
43 minIndex = inda[0][0]
44 except:
45 minIndex = 0
46
47 try:
48 maxIndex = indb[0][-1]
49 except:
50 maxIndex = len(heightList)
51 return minIndex,maxIndex
Show file before | Show file after | Hide comments
@@ -53,18 +53,24 class ProcessingUnit(object):
53 53 return self.operations[objId]
54 54
55 55 def call(self, **kwargs):
56 '''
57 '''
58 56
57 mybool = (self.dataOut.type == 'Voltage') and self.dataOut.useInputBuffer and (not self.dataOut.buffer_empty) #liberar desde buffer
59 58 try:
60 59
61 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
62 return self.dataIn.isReady()
63 elif self.dataIn is None or not self.dataIn.error:
60
61 if mybool:
62 #print("run jeje")
64 63 self.run(**kwargs)
65 elif self.dataIn.error:
66 self.dataOut.error = self.dataIn.error
67 self.dataOut.flagNoData = True
64 else:
65 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
66 return self.dataIn.isReady()
67 elif self.dataIn is None or not self.dataIn.error: #unidad de lectura o procesamiento regular
68 self.run(**kwargs)
69 elif self.dataIn.error:
70 self.dataOut.error = self.dataIn.error
71 self.dataOut.flagNoData = True
72 print("exec proc error")
73
68 74 except:
69 75
70 76 err = traceback.format_exc()
@@ -76,8 +82,10 class ProcessingUnit(object):
76 82 log.error(err, self.name)
77 83 self.dataOut.error = True
78 84
85
79 86 for op, optype, opkwargs in self.operations:
80 if optype == 'other' and self.dataOut.isReady():
87
88 if (optype == 'other' and self.dataOut.isReady()) or mybool:
81 89 try:
82 90 self.dataOut = op.run(self.dataOut, **opkwargs)
83 91 except Exception as e:
@@ -89,7 +97,21 class ProcessingUnit(object):
89 97 elif optype == 'external' and self.dataOut.error:
90 98 op.queue.put(copy.deepcopy(self.dataOut))
91 99
92 return 'Error' if self.dataOut.error else self.dataOut.isReady()
100
101 if not self.dataOut.error:
102 if self.dataOut.type == 'Voltage':
103 if not self.dataOut.buffer_empty : #continue
104 return 'no_Read'
105 elif self.dataOut.useInputBuffer and (self.dataOut.buffer_empty) and self.dataOut.isReady() :
106 return 'new_Read'
107 else:
108 return True
109 else:
110 #print("ret True")
111 return True
112 else:
113 return 'Error'
114 #return 'Error' if self.dataOut.error else True #self.dataOut.isReady()
93 115
94 116 def setup(self):
95 117
Show file before | Show file after | Hide comments
@@ -39,6 +39,8 class SpectraProc(ProcessingUnit):
39 39
40 40 def __updateSpecFromVoltage(self):
41 41
42
43
42 44 self.dataOut.timeZone = self.dataIn.timeZone
43 45 self.dataOut.dstFlag = self.dataIn.dstFlag
44 46 self.dataOut.errorCount = self.dataIn.errorCount
@@ -130,7 +132,7 class SpectraProc(ProcessingUnit):
130 132 self.dataOut.copy(self.dataIn)
131 133
132 134 except Exception as e:
133 print(e)
135 print("Error dataIn ",e)
134 136
135 137 if shift_fft:
136 138 #desplaza a la derecha en el eje 2 determinadas posiciones
@@ -158,7 +160,7 class SpectraProc(ProcessingUnit):
158 160 self.dataOut.ippFactor = 1
159 161
160 162 self.dataOut.nFFTPoints = nFFTPoints
161
163 #print(" volts ch,prof, h: ", self.dataIn.data.shape)
162 164 if self.buffer is None:
163 165 self.buffer = numpy.zeros((self.dataIn.nChannels,
164 166 nProfiles,
@@ -195,6 +197,7 class SpectraProc(ProcessingUnit):
195 197 self.firstdatatime = self.dataIn.utctime
196 198
197 199 if self.profIndex == nProfiles:
200
198 201 self.__updateSpecFromVoltage()
199 202 if pairsList == None:
200 203 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
@@ -209,6 +212,7 class SpectraProc(ProcessingUnit):
209 212 raise ValueError("The type of input object '%s' is not valid".format(
210 213 self.dataIn.type))
211 214
215
212 216 def __selectPairs(self, pairsList):
213 217
214 218 if not pairsList:
@@ -492,14 +496,14 class removeDC(Operation):
492 496 return self.dataOut
493 497
494 498 class getNoise(Operation):
495
499 warnings = False
496 500 def __init__(self):
497 501
498 502 Operation.__init__(self)
499 503
500 def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None):
504 def run(self, dataOut, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
501 505 self.dataOut = dataOut
502
506 self.warnings = warnings
503 507 if minHei == None:
504 508 minHei = self.dataOut.heightList[0]
505 509
@@ -507,13 +511,15 class getNoise(Operation):
507 511 maxHei = self.dataOut.heightList[-1]
508 512
509 513 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
510 print('minHei: %.2f is out of the heights range' % (minHei))
511 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
514 if self.warnings:
515 print('minHei: %.2f is out of the heights range' % (minHei))
516 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
512 517 minHei = self.dataOut.heightList[0]
513 518
514 519 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
515 print('maxHei: %.2f is out of the heights range' % (maxHei))
516 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
520 if self.warnings:
521 print('maxHei: %.2f is out of the heights range' % (maxHei))
522 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
517 523 maxHei = self.dataOut.heightList[-1]
518 524
519 525
@@ -535,13 +541,15 class getNoise(Operation):
535 541 maxFreq = freqrange[-1]
536 542
537 543 if (minFreq < freqrange[0]) or (minFreq > maxFreq):
538 print('minFreq: %.2f is out of the frequency range' % (minFreq))
539 print('minFreq is setting to %.2f' % (freqrange[0]))
544 if self.warnings:
545 print('minFreq: %.2f is out of the frequency range' % (minFreq))
546 print('minFreq is setting to %.2f' % (freqrange[0]))
540 547 minFreq = freqrange[0]
541 548
542 549 if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
543 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
544 print('maxFreq is setting to %.2f' % (freqrange[-1]))
550 if self.warnings:
551 print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
552 print('maxFreq is setting to %.2f' % (freqrange[-1]))
545 553 maxFreq = freqrange[-1]
546 554
547 555 indminPoint = numpy.where(freqrange >= minFreq)
@@ -557,13 +565,15 class getNoise(Operation):
557 565 maxVel = velrange[-1]
558 566
559 567 if (minVel < velrange[0]) or (minVel > maxVel):
560 print('minVel: %.2f is out of the velocity range' % (minVel))
561 print('minVel is setting to %.2f' % (velrange[0]))
568 if self.warnings:
569 print('minVel: %.2f is out of the velocity range' % (minVel))
570 print('minVel is setting to %.2f' % (velrange[0]))
562 571 minVel = velrange[0]
563 572
564 573 if (maxVel > velrange[-1]) or (maxVel < minVel):
565 print('maxVel: %.2f is out of the velocity range' % (maxVel))
566 print('maxVel is setting to %.2f' % (velrange[-1]))
574 if self.warnings:
575 print('maxVel: %.2f is out of the velocity range' % (maxVel))
576 print('maxVel is setting to %.2f' % (velrange[-1]))
567 577 maxVel = velrange[-1]
568 578
569 579 indminPoint = numpy.where(velrange >= minVel)
@@ -1094,14 +1104,14 class IntegrationFaradaySpectra(Operation):
1094 1104 n = None
1095 1105 minHei_ind = None
1096 1106 maxHei_ind = None
1097 avg = 1.0
1107 navg = 1.0
1098 1108 factor = 0.0
1099 1109
1100 1110 def __init__(self):
1101 1111
1102 1112 Operation.__init__(self)
1103 1113
1104 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, factor=0.75):
1114 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
1105 1115 """
1106 1116 Set the parameters of the integration class.
1107 1117
@@ -1125,6 +1135,7 class IntegrationFaradaySpectra(Operation):
1125 1135 self.__byTime = False
1126 1136
1127 1137 self.factor = factor
1138 self.navg = avg
1128 1139 #self.ByLags = dataOut.ByLags ###REDEFINIR
1129 1140 self.ByLags = False
1130 1141
@@ -1277,13 +1288,11 class IntegrationFaradaySpectra(Operation):
1277 1288 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1278 1289 indexes=numpy.array(indexes)
1279 1290 indexmin=numpy.min(indexes)
1280 #print("clean CH: ", i)
1291
1281 1292 if indexmin != buffer1.shape[0]:
1282 1293 if self.nChannels > 1:
1283 1294 cspc_outliers_exist= True
1284 #print("outliers cspc")
1285 ###sortdata=numpy.sort(buffer1,axis=0)
1286 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1295
1287 1296 lt=outliers_IDs
1288 1297 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1289 1298
@@ -1400,7 +1409,7 class IntegrationFaradaySpectra(Operation):
1400 1409
1401 1410 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1402 1411
1403 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, factor=0.75):
1412 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
1404 1413 self.dataOut = dataOut.copy()
1405 1414 if n == 1:
1406 1415 return self.dataOut
@@ -1410,7 +1419,7 class IntegrationFaradaySpectra(Operation):
1410 1419 self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
1411 1420 #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc)
1412 1421 if not self.isConfig:
1413 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, factor)
1422 self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
1414 1423 self.isConfig = True
1415 1424
1416 1425 if not self.ByLags:
@@ -1742,7 +1751,9 class IncohInt(Operation):
1742 1751 Add a profile to the __buffer_spc and increase in one the __profileIndex
1743 1752
1744 1753 """
1745
1754 if data_spc.all() == numpy.nan :
1755 print("nan ")
1756 return
1746 1757 self.__buffer_spc += data_spc
1747 1758
1748 1759 if data_cspc is None:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (524 lines changed) Show them Hide them
@@ -4,9 +4,11 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 6 from schainpy.utils import log
7 from schainpy.model.io.utils import getHei_index
7 8 from time import time
9 import datetime
8 10 import numpy
9
11 import copy
10 12
11 13 class VoltageProc(ProcessingUnit):
12 14
@@ -19,12 +21,15 class VoltageProc(ProcessingUnit):
19 21 self.setupReq = False
20 22
21 23 def run(self):
22
24 #print("running volt proc")
23 25 if self.dataIn.type == 'AMISR':
24 26 self.__updateObjFromAmisrInput()
25 27
26 if self.dataIn.type == 'Voltage':
27 self.dataOut.copy(self.dataIn)
28 if self.dataOut.buffer_empty:
29 if self.dataIn.type == 'Voltage':
30 self.dataOut.copy(self.dataIn)
31 #print("new volts reading")
32
28 33
29 34 def __updateObjFromAmisrInput(self):
30 35
@@ -1680,6 +1685,7 class SSheightProfiles(Operation):
1680 1685 dataOut.flagNoData = True
1681 1686
1682 1687 profileIndex = None
1688 #print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights)
1683 1689 #print(dataOut.getFreqRange(1)/1000.)
1684 1690 #exit(1)
1685 1691 if dataOut.flagDataAsBlock:
@@ -1691,30 +1697,11 class SSheightProfiles(Operation):
1691 1697 #print("Setup done")
1692 1698 self.isConfig = True
1693 1699
1694 #DC_Hae = numpy.array([0.398+0.588j, -0.926+0.306j, -0.536-0.682j, -0.072+0.53j, 0.368-0.356j, 0.996+0.362j])
1695 #DC_Hae = numpy.array([ 0.001025 +0.0516375j, 0.03485 +0.20923125j, -0.168 -0.02720625j,
1696 #-0.1105375 +0.0707125j, -0.20309375-0.09670625j, 0.189775 +0.02716875j])*(-3.5)
1697
1698 #DC_Hae = numpy.array([ -32.26 +8.66j, -32.26 +8.66j])
1699 1700
1700 #DC_Hae = numpy.array([-2.78500000e-01 -1.39175j, -6.63237294e+02+210.4268625j])
1701
1702 #print(dataOut.data[0,13:15])
1703 #dataOut.data = dataOut.data - DC_Hae[:,None]
1704 #print(dataOut.data[0,13:15])
1705 #exit(1)
1706 1701 if code is not None:
1707 1702 code = numpy.array(code)
1708 1703 code_block = code
1709 '''
1710 roll = 0
1711 code = numpy.roll(code,roll,axis=0)
1712 code = numpy.reshape(code,(5,100,64))
1713 block = dataOut.CurrentBlock%5
1714
1715 day_dif = 0 #day_19_Oct_2021: 3
1716 code_block = code[block-1-day_dif,:,:]
1717 '''
1704
1718 1705 if repeat is not None:
1719 1706 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1720 1707 #print(code_block.shape)
@@ -1756,3 +1743,492 class SSheightProfiles(Operation):
1756 1743 #exit(1)
1757 1744
1758 1745 return dataOut
1746 ################################################################################3############################3
1747 ################################################################################3############################3
1748 ################################################################################3############################3
1749 ################################################################################3############################3
1750
1751 class SSheightProfiles2(Operation):
1752 '''
1753 Procesa por perfiles y por bloques
1754 '''
1755
1756 step = None
1757 nsamples = None
1758 bufferShape = None
1759 profileShape = None
1760 sshProfiles = None
1761 profileIndex = None
1762
1763 def __init__(self, **kwargs):
1764
1765 Operation.__init__(self, **kwargs)
1766 self.isConfig = False
1767
1768 def setup(self,dataOut ,step = None , nsamples = None):
1769
1770 if step == None and nsamples == None:
1771 raise ValueError("step or nheights should be specified ...")
1772
1773 self.step = step
1774 self.nsamples = nsamples
1775 self.__nChannels = int(dataOut.nChannels)
1776 self.__nProfiles = int(dataOut.nProfiles)
1777 self.__nHeis = int(dataOut.nHeights)
1778
1779 residue = (self.__nHeis - self.nsamples) % self.step
1780 if residue != 0:
1781 print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
1782
1783 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1784 #numberProfile = self.nsamples
1785 self.numberSamples = (self.__nHeis - self.nsamples)/self.step
1786
1787 self.new_nHeights = self.numberSamples
1788
1789 self.bufferShape = int(self.__nChannels), int(self.numberSamples), int(self.nsamples) # nchannels, nsamples , nprofiles
1790 self.profileShape = int(self.__nChannels), int(self.nsamples), int(self.numberSamples) # nchannels, nprofiles, nsamples
1791
1792 self.buffer = numpy.zeros(self.bufferShape , dtype=numpy.complex)
1793 self.sshProfiles = numpy.zeros(self.profileShape, dtype=numpy.complex)
1794
1795 def getNewProfiles(self, data, code=None, repeat=None):
1796
1797 if code is not None:
1798 code = numpy.array(code)
1799 code_block = code
1800
1801 if repeat is not None:
1802 code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
1803 #print("buff, data, :",self.buffer.shape, data.shape)
1804 for i in range(int(self.new_nHeights)): #nuevas alturas
1805 if code is not None:
1806 self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]*code_block
1807 else:
1808 self.buffer[:,i] = data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
1809
1810 for j in range(self.__nChannels): #en los cananles
1811 self.sshProfiles[j] = numpy.transpose(self.buffer[j])
1812
1813
1814
1815
1816 def run(self, dataOut, step, nsamples, code = None, repeat = None):
1817
1818 dataOut.flagNoData = True
1819 #print("init data shape:", dataOut.data.shape)
1820 #print("ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
1821 # int(dataOut.nProfiles),int(dataOut.nHeights)))
1822
1823 profileIndex = None
1824 # if not dataOut.flagDataAsBlock:
1825 # dataOut.nProfiles = 1
1826
1827 if not self.isConfig:
1828 self.setup(dataOut, step=step , nsamples=nsamples)
1829 #print("Setup done")
1830 self.isConfig = True
1831
1832 dataBlock = None
1833
1834 nprof = 1
1835 if dataOut.flagDataAsBlock:
1836 nprof = int(dataOut.nProfiles)
1837
1838 #print("dataOut nProfiles:", dataOut.nProfiles)
1839 for profile in range(nprof):
1840 if dataOut.flagDataAsBlock:
1841 self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
1842 else:
1843 self.getNewProfiles(dataOut.data, code=code, repeat=repeat)
1844 if profile == 0:
1845 dataBlock = self.sshProfiles.copy()
1846 else: #by blocks
1847 dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
1848 print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
1849
1850 profileIndex = self.nsamples
1851 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1852 ippSeconds = (deltaHeight*1.0e-6)/(0.15)
1853
1854
1855 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)
1859 dataOut.ippSeconds = ippSeconds
1860 dataOut.step = self.step
1861 dataOut.flagNoData = False
1862 if dataOut.flagDataAsBlock:
1863 dataOut.nProfiles = int(dataOut.nProfiles*self.nsamples)
1864
1865 else:
1866 dataOut.nProfiles = int(self.nsamples)
1867 dataOut.profileIndex = dataOut.nProfiles
1868 dataOut.flagDataAsBlock = True
1869
1870 dataBlock = None
1871 #print("new data shape:", dataOut.data.shape)
1872
1873 return dataOut
1874
1875
1876
1877
1878 #import skimage.color
1879 #import skimage.io
1880 import matplotlib.pyplot as plt
1881
1882 class clean2DArray(Operation):
1883 '''
1884
1885 '''
1886 isConfig = False
1887 n = None
1888 __timeInterval = None
1889 __profIndex = 0
1890 __byTime = False
1891 __dataReady = False
1892 __buffer_data = []
1893 __buffer_times = []
1894 __initime = None
1895 __count_exec = 0
1896 __profIndex = 0
1897 buffer = None
1898 lenProfileOut = 1
1899
1900 init_prof = 0
1901 end_prof = 0
1902 n_profiles = 0
1903 first_utcBlock = None
1904 __dh = 0
1905
1906 def __init__(self, **kwargs):
1907
1908 Operation.__init__(self, **kwargs)
1909 self.isConfig = False
1910
1911 def setup(self,dataOut, n=None , timeInterval=None):
1912
1913 if n == None and timeInterval == None:
1914 raise ValueError("nprofiles or timeInterval should be specified ...")
1915
1916 if n != None:
1917 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
1924 self.__profIndex = 0
1925 self.buffer = None
1926 self.lenProfileOut = 1
1927
1928 self.init_prof = 0
1929 self.end_prof = 0
1930 self.n_profiles = 0
1931 self.first_utcBlock = None
1932 self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
1933
1934 def cleanOutliersByBlock(self):
1935 #print(self.__buffer_data[0].shape)
1936 data = self.__buffer_data#.copy()
1937 #print("cleaning shape inpt: ",data.shape)
1938
1939 self.__buffer_data = []
1940
1941 spectrum = numpy.fft.fft2(data, axes=(0,2))
1942 #print("spc : ",spectrum.shape)
1943 (nch,nsamples, nh) = spectrum.shape
1944
1945 #nch =
1946
1947 for ch in range(nch):
1948 dh = self.__dh
1949 dt1 = (dh*1.0e-6)/(0.15)
1950 dt2 = self.__buffer_times[1]-self.__buffer_times[0]
1951 freqv = numpy.fft.fftfreq(nh, d=dt1)
1952 freqh = numpy.fft.fftfreq(self.n, d=dt2)
1953 #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
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
1971 plt.show()
1972
1973
1974
1975 #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
1976
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
1998
1999 if self.__profIndex == self.n:
2000 #print("apnd : ",data)
2001 dataBlock = self.cleanOutliersByBlock()
2002 self.__dataReady = True
2003
2004 return dataBlock
2005
2006 def fillBuffer(self, data, datatime):
2007
2008 if self.__profIndex == 0:
2009 self.__buffer_data = data.copy()
2010
2011 else:
2012 self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
2013 self.__profIndex += 1
2014 self.__buffer_times.append(datatime)
2015
2016 def getData(self, data, datatime=None):
2017
2018 if self.__profIndex == 0:
2019 self.__initime = datatime
2020
2021 if self.__byTime:
2022 dataBlock = self.byTime(data, datatime)
2023 else:
2024 dataBlock = self.byProfiles(data, datatime)
2025
2026
2027
2028 if dataBlock is None:
2029 return None
2030
2031
2032
2033 return dataBlock
2034
2035 def releaseBlock(self, dataOut):
2036
2037 if self.n % self.lenProfileOut != 0:
2038 raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n_profiles))
2039 return None
2040
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
2044 self.init_prof = self.end_prof
2045 self.end_prof += self.lenProfileOut
2046 print("data release shape: ",dataOut.data.shape, self.end_prof)
2047
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
2055 def run(self, dataOut, n=None, timeInterval=None, nProfilesOut=1):
2056 #print("run op buffer 2D")
2057 self.nChannels = dataOut.nChannels
2058 self.nHeights = dataOut.nHeights
2059
2060 if not self.isConfig:
2061 self.setup(dataOut,n=n, timeInterval=timeInterval)
2062 self.isConfig = True
2063
2064 dataBlock = None
2065
2066 if not dataOut.buffer_empty: #hay datos acumulados
2067
2068 if self.init_prof == 0:
2069 self.lenProfileOut = nProfilesOut
2070 dataOut.flagNoData = False
2071 #print("tp 2 ",dataOut.data.shape)
2072 #(ch, self.n_profiles, nh) = self.buffer.shape
2073 #print("tp 3 ",self.dataOut.data.shape)
2074 #print("rel: ",self.buffer[:,-1,:])
2075 self.init_prof = 0
2076 self.end_prof = self.lenProfileOut
2077 self.first_utcBlock = dataOut.utctime
2078 dataOut.nProfiles = self.lenProfileOut
2079 dataOut.error = False
2080 dataOut.flagDataAsBlock = True
2081 #print("prof: ",self.init_prof)
2082 dataOut.flagNoData = False
2083 return self.releaseBlock(dataOut)
2084
2085
2086 #print("tp 223 ",dataOut.data.shape)
2087 dataOut.flagNoData = True
2088
2089
2090
2091 try:
2092 #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
2093 dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
2094 self.__count_exec +=1
2095 except Exception as e:
2096 print("Error getting profiles data",self.__count_exec )
2097 print(e)
2098 sys.exit()
2099
2100 if self.__dataReady:
2101 self.__count_exec = 0
2102 #dataOut.data =
2103 self.buffer = numpy.flip(dataBlock, axis=1)
2104
2105 dataOut.utctime = self.__initime
2106 dataOut.nProfiles = self.__profIndex
2107 #dataOut.flagNoData = False
2108 self.__profIndex = 0
2109 self.__initime = None
2110 dataBlock = None
2111 self.__buffer_times = []
2112 dataOut.error = False
2113 dataOut.useInputBuffer = True
2114 dataOut.buffer_empty = False
2115 print("1 ch: {} prof: {} hs: {}".format(int(dataOut.nChannels),
2116 int(dataOut.nProfiles),int(dataOut.nHeights)))
2117 #return None
2118
2119
2120 #print(self.__count_exec)
2121
2122 return dataOut
2123
2124 class CleanProfileSats(Operation):
2125 '''
2126 Omite los perfiles contaminados con señal de satelites,
2127 In: minHei = min_sat_range
2128 max_sat_range
2129 min_hei_ref
2130 max_hei_ref
2131 th = diference between profiles mean, ref and sats
2132 Out:
2133 profile clean
2134 '''
2135
2136 isConfig = False
2137 min_sats = 0
2138 max_sats = 999999999
2139 min_ref= 0
2140 max_ref= 9999999999
2141 needReshape = False
2142 count = 0
2143 thdB = 0
2144 byRanges = False
2145 min_sats = None
2146 max_sats = None
2147
2148 def __init__(self, **kwargs):
2149
2150 Operation.__init__(self, **kwargs)
2151 self.isConfig = False
2152
2153
2154 def setup(self, dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList):
2155
2156 if rangeHeiList!=None:
2157 self.byRanges = True
2158 else:
2159 if minHei==None or maxHei==None :
2160 raise ValueError("Parameters heights are required")
2161 if minRef==None or maxRef==None:
2162 raise ValueError("Parameters heights are required")
2163
2164 if self.byRanges:
2165 self.min_sats = []
2166 self.max_sats = []
2167 for min,max in rangeHeiList:
2168 a,b = getHei_index(min, max, dataOut.heightList)
2169 self.min_sats.append(a)
2170 self.max_sats.append(b)
2171 else:
2172 self.min_sats, self.max_sats = getHei_index(minHei, maxHei, dataOut.heightList)
2173 self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
2174 self.th = th
2175 self.thdB = thdB
2176 self.isConfig = True
2177
2178
2179 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)
2183
2184 sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
2185 p_sats = 10*numpy.log10(sats.real)
2186 m_sats = numpy.mean(p_sats)
2187
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))
2190 #print("Removing profiles...")
2191 return False
2192
2193 return True
2194
2195 def isProfileClean(self, data):
2196 '''
2197 Analiza solo 1 canal, y descarta todos...
2198 '''
2199
2200 clean = True
2201
2202 if self.byRanges:
2203
2204 for n in range(len(self.min_sats)):
2205 c = self.compareRanges(data,self.min_sats[n],self.max_sats[n])
2206 clean = clean and c
2207 else:
2208
2209 clean = (self.compareRanges(data, self.min_sats,self.max_sats))
2210
2211 return clean
2212
2213
2214
2215 def run(self, dataOut, minHei=None, maxHei=None, minRef=None, maxRef=None, th=5, thdB=65, rangeHeiList=None):
2216 dataOut.flagNoData = True
2217
2218 if not self.isConfig:
2219 self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
2220 self.isConfig = True
2221
2222 if dataOut.flagDataAsBlock:
2223 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
2224
2225 else:
2226
2227 if not self.isProfileClean(dataOut.data):
2228 return dataOut
2229 #dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)
2230 #self.count += 1
2231
2232 dataOut.flagNoData = False
2233
2234 return dataOut
General Comments 0
You need to be logged in to leave comments. Login now