##// END OF EJS Templates
schain-jc
RSS
Review last commit
Juan C. Espinoza -
r1186:29f68a738921 merge
parent child
Show More
Side by Side Unified
Context file:
r1186:29f68a738921 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,200 +1,208
1 1 import click
2 2 import schainpy
3 3 import subprocess
4 4 import os
5 5 import sys
6 6 import glob
7 7 save_stdout = sys.stdout
8 8 sys.stdout = open('/dev/null', 'w')
9 9 from multiprocessing import cpu_count
10 10 from schainpy.controller import Project
11 11 from schainpy.model import Operation, ProcessingUnit
12 12 from schainpy.utils import log
13 13 from importlib import import_module
14 14 from pydoc import locate
15 15 from fuzzywuzzy import process
16 16 from schainpy.cli import templates
17 17 sys.stdout = save_stdout
18 18
19 19
20 20 def getProcs():
21 21 modules = dir(schainpy.model)
22 22 procs = check_module(modules, ProcessingUnit)
23 23 try:
24 24 procs.remove('ProcessingUnit')
25 25 except Exception as e:
26 26 pass
27 27 return procs
28 28
29 29 def getOperations():
30 30 module = dir(schainpy.model)
31 31 noProcs = [x for x in module if not x.endswith('Proc')]
32 32 operations = check_module(noProcs, Operation)
33 33 try:
34 34 operations.remove('Operation')
35 35 except Exception as e:
36 36 pass
37 37 return operations
38 38
39 39 def getArgs(op):
40 40 module = locate('schainpy.model.{}'.format(op))
41 41 args = module().getAllowedArgs()
42 42 try:
43 43 args.remove('self')
44 44 except Exception as e:
45 45 pass
46 46 try:
47 47 args.remove('dataOut')
48 48 except Exception as e:
49 49 pass
50 50 return args
51 51
52 52 def getAll():
53 53 allModules = dir(schainpy.model)
54 54 modules = check_module(allModules, Operation)
55 55 modules.extend(check_module(allModules, ProcessingUnit))
56 56 return modules
57 57
58 58
59 59 def print_version(ctx, param, value):
60 60 if not value or ctx.resilient_parsing:
61 61 return
62 62 click.echo(schainpy.__version__)
63 63 ctx.exit()
64 64
65 65
66 66 PREFIX = 'experiment'
67 67
68 68 @click.command()
69 69 @click.option('--version', '-v', is_flag=True, callback=print_version, help='SChain version', type=str)
70 70 @click.argument('command', default='run', required=True)
71 71 @click.argument('nextcommand', default=None, required=False, type=str)
72 72 def main(command, nextcommand, version):
73 """COMMAND LINE INTERFACE FOR SIGNAL CHAIN - JICAMARCA RADIO OBSERVATORY \n
73 """COMMAND LINE INTERFACE FOR SIGNAL CHAIN - JICAMARCA RADIO OBSERVATORY V3.0\n
74 74 Available commands.\n
75 --xml: runs a schain XML generated file\n
75 xml: runs a schain XML generated file\n
76 76 run: runs any python script starting 'experiment_'\n
77 77 generate: generates a template schain script\n
78 78 search: return avilable operations, procs or arguments of the give operation/proc\n"""
79 79 if command == 'xml':
80 80 runFromXML(nextcommand)
81 81 elif command == 'generate':
82 82 generate()
83 83 elif command == 'test':
84 84 test()
85 85 elif command == 'run':
86 86 runschain(nextcommand)
87 87 elif command == 'search':
88 88 search(nextcommand)
89 89 else:
90 90 log.error('Command {} is not defined'.format(command))
91 91
92 92
93 93 def check_module(possible, instance):
94 94 def check(x):
95 95 try:
96 96 instancia = locate('schainpy.model.{}'.format(x))
97 97 return isinstance(instancia(), instance)
98 98 except Exception as e:
99 99 return False
100 100 clean = clean_modules(possible)
101 101 return [x for x in clean if check(x)]
102 102
103 103
104 104 def clean_modules(module):
105 105 noEndsUnder = [x for x in module if not x.endswith('__')]
106 106 noStartUnder = [x for x in noEndsUnder if not x.startswith('__')]
107 107 noFullUpper = [x for x in noStartUnder if not x.isupper()]
108 108 return noFullUpper
109 109
110 110
111 111 def search(nextcommand):
112 112 if nextcommand is None:
113 113 log.error('There is no Operation/ProcessingUnit to search', '')
114 114 elif nextcommand == 'procs':
115 115 procs = getProcs()
116 116 log.success(
117 117 'Current ProcessingUnits are:\n{}'.format('\n'.join(procs)), '')
118 118
119 119 elif nextcommand == 'operations':
120 120 operations = getOperations()
121 121 log.success('Current Operations are:\n{}'.format(
122 122 '\n'.join(operations)), '')
123 123 else:
124 124 try:
125 125 args = getArgs(nextcommand)
126 126 if len(args) == 0:
127 127 log.success('`{}` has no arguments'.format(nextcommand), '')
128 128 else:
129 129 log.success('`{}` arguments: {}'.format(
130 130 nextcommand, ', '.join(args)), '')
131 131 except Exception as e:
132 132 log.error('Module `{}` does not exists'.format(nextcommand), '')
133 133 allModules = getAll()
134 134 similar = [t[0] for t in process.extract(nextcommand, allModules, limit=12) if t[1]>80]
135 135 log.success('Possible modules are: {}'.format(', '.join(similar)), '')
136 136
137 137 def runschain(nextcommand):
138 138 if nextcommand is None:
139 139 currentfiles = glob.glob('./{}_*.py'.format(PREFIX))
140 140 numberfiles = len(currentfiles)
141 141 if numberfiles > 1:
142 142 log.error('There is more than one file to run')
143 143 elif numberfiles == 1:
144 144 subprocess.call(['python ' + currentfiles[0]], shell=True)
145 145 else:
146 146 log.error('There is no file to run')
147 147 else:
148 148 try:
149 149 subprocess.call(['python ' + nextcommand], shell=True)
150 150 except Exception as e:
151 151 log.error("I cannot run the file. Does it exists?")
152 152
153 153
154 154 def basicInputs():
155 155 inputs = {}
156 inputs['desc'] = click.prompt(
157 'Enter a description', default="A schain project", type=str)
158 156 inputs['name'] = click.prompt(
159 157 'Name of the project', default="project", type=str)
158 inputs['desc'] = click.prompt(
159 'Enter a description', default="A schain project", type=str)
160 inputs['multiprocess'] = click.prompt(
161 '''Select data type:
162
163 - Voltage (*.r): [1]
164 - Spectra (*.pdata): [2]
165 - Voltage and Spectra (*.r): [3]
166
167 -->''', type=int)
160 168 inputs['path'] = click.prompt('Data path', default=os.getcwd(
161 169 ), type=click.Path(exists=True, resolve_path=True))
162 170 inputs['startDate'] = click.prompt(
163 171 'Start date', default='1970/01/01', type=str)
164 172 inputs['endDate'] = click.prompt(
165 'End date', default='2017/12/31', type=str)
173 'End date', default='2018/12/31', type=str)
166 174 inputs['startHour'] = click.prompt(
167 175 'Start hour', default='00:00:00', type=str)
168 176 inputs['endHour'] = click.prompt('End hour', default='23:59:59', type=str)
169 177 inputs['figpath'] = inputs['path'] + '/figs'
170 178 return inputs
171 179
172 180
173 181 def generate():
174 182 inputs = basicInputs()
175 inputs['multiprocess'] = click.confirm('Is this a multiprocess script?')
176 if inputs['multiprocess']:
177 inputs['nProcess'] = click.prompt(
178 'How many process?', default=cpu_count(), type=int)
179 current = templates.multiprocess.format(**inputs)
180 else:
181 current = templates.basic.format(**inputs)
183
184 if inputs['multiprocess'] == 1:
185 current = templates.voltage.format(**inputs)
186 elif inputs['multiprocess'] == 2:
187 current = templates.spectra.format(**inputs)
188 elif inputs['multiprocess'] == 3:
189 current = templates.voltagespectra.format(**inputs)
182 190 scriptname = '{}_{}.py'.format(PREFIX, inputs['name'])
183 191 script = open(scriptname, 'w')
184 192 try:
185 193 script.write(current)
186 194 log.success('Script {} generated'.format(scriptname))
187 195 except Exception as e:
188 196 log.error('I cannot create the file. Do you have writing permissions?')
189 197
190 198
191 199 def test():
192 200 log.warning('testing')
193 201
194 202
195 203 def runFromXML(filename):
196 204 controller = Project()
197 205 if not controller.readXml(filename):
198 206 return
199 207 controller.start()
200 208 return
Show file before | Show file after | Hide comments
@@ -1,90 +1,269
1 basic = '''from schainpy.controller import Project
1 voltage = '''import os, sys, time
2 from schainpy.controller import Project
3
4
5 def main():
6 desc = "{desc}"
7 controller = Project()
8 controller.setup(id='200', name="{name}", description=desc)
9
10 read_unit = controller.addReadUnit(datatype='Voltage',
11 path="{path}",
12 startDate="{startDate}",
13 endDate="{endDate}",
14 startTime="{startHour}",
15 endTime="{endHour}",
16 online=0,
17 verbose=1,
18 walk=0,
19 delay=180,
20 )
21
22 code = '[[1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1]]'
23 nCode = '128'
24 nBaud = '3'
25
26
27 proc_voltage = controller.addProcUnit(name='VoltageProc', inputId=read_unit.getId())
28
29 op1 = proc_voltage.addOperation(name='selectChannels', optype='self')
30 op1.addParameter(name='channelList', value='0, 1, 2, 3', format='intlist')
31
32 op2 = proc_voltage.addOperation(name='filterByHeights', optype='self')
33 op2.addParameter(name='window', value='4', format='int')
34
35 op3 = proc_voltage.addOperation(name='ProfileSelector', optype='other')
36 op3.addParameter(name='profileRangeList', value='32, 159', format='intList')
37
38 op4 = proc_voltage.addOperation(name='Decoder', optype='other')
39 op4.addParameter(name='code', value=code, format='intlist')
40 op4.addParameter(name='nCode', value=nCode, format='int')
41 op4.addParameter(name='nBaud', value=nBaud, format='int')
42 op4.addParameter(name='mode', value='0', format='int')
43
44 op5 = proc_voltage.addOperation(name='Scope', optype='external')
45 op5.addParameter(name='id', value='30', format='int')
46
47
48
49
50
51 controller.start()
52
53 if __name__ == '__main__':
54 import time
55 start_time = time.time()
56 main()
57 print("--- %s seconds ---" % (time.time() - start_time))
2 58
3 desc = "{desc}"
4 project = Project()
5 project.setup(id='200', name="{name}", description=desc)
6
7 voltage_reader = project.addReadUnit(datatype='VoltageReader',
8 path="{path}",
9 startDate="{startDate}",
10 endDate="{endDate}",
11 startTime="{startHour}",
12 endTime="{endHour}",
13 online=0,
14 verbose=1,
15 walk=1,
16 )
17
18 voltage_proc = project.addProcUnit(datatype='VoltageProc', inputId=voltage_reader.getId())
19
20 profile = voltage_proc.addOperation(name='ProfileSelector', optype='other')
21 profile.addParameter(name='profileRangeList', value='120,183', format='intlist')
22
23 rti = voltage_proc.addOperation(name='RTIPlot', optype='other')
24 rti.addParameter(name='wintitle', value='Jicamarca Radio Observatory', format='str')
25 rti.addParameter(name='showprofile', value='0', format='int')
26 rti.addParameter(name='xmin', value='0', format='int')
27 rti.addParameter(name='xmax', value='24', format='int')
28 rti.addParameter(name='figpath', value="{figpath}", format='str')
29 rti.addParameter(name='wr_period', value='5', format='int')
30 rti.addParameter(name='exp_code', value='22', format='int')
31
32
33 project.start()
34 59 '''
35 60
61
62 spectra = '''import os, sys, time
63 from schainpy.controller import Project
64
65
66 def main():
67 desc = "{desc}"
68 controller = Project()
69 controller.setup(id='300', name="{name}", description=desc)
70
71 read_unit = controller.addReadUnit(datatype='Spectra',
72 path="{path}",
73 startDate="{startDate}",
74 endDate="{endDate}",
75 startTime="{startHour}",
76 endTime="{endHour}",
77 online=0,
78 verbose=1,
79 walk=0,
80 delay=180,
81 )
82
83 proc_spectra = controller.addProcUnit(datatype='Spectra', inputId=read_unit.getId())
84 proc_spectra.addParameter(name='nFFTPoints', value='128', format='int')
85 proc_spectra.addParameter(name='nProfiles', value='128', format='int')
86 proc_spectra.addParameter(name='pairsList', value='(0, 1), (2, 3)', format='pairslist')
87
88 op1 = proc_spectra.addOperation(name='IncohInt', optype='other')
89 op1.addParameter(name='n', value='4', format='int')
90
91 op2 = proc_spectra.addOperation(name='CrossSpectraPlot', optype='external')
92 op2.addParameter(name='id', value='10', format='int')
93 op2.addParameter(name='zmin', value='10.0', format='float')
94 op2.addParameter(name='zmax', value='35.0', format='float')
95
96
97 op3 = proc_spectra.addOperation(name='RTIPlot', optype='external')
98 op3.addParameter(name='id', value='20', format='int')
99 op3.addParameter(name='wintitle', value='RTI', format='str')
100 op3.addParameter(name='xmin', value='0', format='float')
101 op3.addParameter(name='xmax', value='24', format='float')
102 op3.addParameter(name='zmin', value='12', format='int')
103 op3.addParameter(name='zmax', value='32', format='int')
104 op3.addParameter(name='showprofile', value='1', format='int')
105 op3.addParameter(name='timerange', value=str(24*60*60), format='int')
106
107 op4 = proc_spectra.addOperation(name='CoherenceMap', optype='external')
108 op4.addParameter(name='id', value='30', format='int')
109 op4.addParameter(name='xmin', value='0.0', format='float')
110 op4.addParameter(name='xmax', value='24.0', format='float')
111
112
113 controller.start()
114
115 if __name__ == '__main__':
116 import time
117 start_time = time.time()
118 main()
119 print("--- %s seconds ---" % (time.time() - start_time))
120
121 '''
122
123 voltagespectra = '''import os, sys, time
124 from schainpy.controller import Project
125
126
127 def main():
128 desc = "{desc}"
129 controller = Project()
130 controller.setup(id='400', name="{name}", description=desc)
131
132 read_unit = controller.addReadUnit(datatype='Voltage',
133 path="{path}",
134 startDate="{startDate}",
135 endDate="{endDate}",
136 startTime="{startHour}",
137 endTime="{endHour}",
138 online=0,
139 verbose=1,
140 walk=0,
141 delay=180,
142 )
143
144 code = '[[1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [1, 1, -1], [-1, -1, 1], [-1, -1, 1], [-1, -1, 1], [1, 1, -1]]'
145 nCode = '128'
146 nBaud = '3'
147
148
149 proc_voltage = controller.addProcUnit(name='VoltageProc', inputId=read_unit.getId())
150
151 op1 = proc_voltage.addOperation(name='selectChannels', optype='self')
152 op1.addParameter(name='channelList', value='0, 1, 2, 3', format='intlist')
153
154 op2 = proc_voltage.addOperation(name='filterByHeights', optype='self')
155 op2.addParameter(name='window', value='4', format='int')
156
157 op3 = proc_voltage.addOperation(name='ProfileSelector', optype='other')
158 op3.addParameter(name='profileRangeList', value='32, 159', format='intList')
159
160 op4 = proc_voltage.addOperation(name='Decoder', optype='other')
161 op4.addParameter(name='code', value=code, format='intlist')
162 op4.addParameter(name='nCode', value=nCode, format='int')
163 op4.addParameter(name='nBaud', value=nBaud, format='int')
164 op4.addParameter(name='mode', value='0', format='int')
165
166
167
168 proc_spectra = controller.addProcUnit(datatype='Spectra', inputId=proc_voltage.getId())
169 proc_spectra.addParameter(name='nFFTPoints', value='128', format='int')
170 proc_spectra.addParameter(name='nProfiles', value='128', format='int')
171 proc_spectra.addParameter(name='pairsList', value='(0, 1), (2, 3)', format='pairslist')
172
173 op5 = proc_spectra.addOperation(name='IncohInt', optype='other')
174 op5.addParameter(name='n', value='4', format='int')
175
176 op6 = proc_spectra.addOperation(name='CrossSpectraPlot', optype='external')
177 op6.addParameter(name='id', value='10', format='int')
178 op6.addParameter(name='zmin', value='10.0', format='float')
179 op6.addParameter(name='zmax', value='35.0', format='float')
180
181
182 op7 = proc_spectra.addOperation(name='RTIPlot', optype='external')
183 op7.addParameter(name='id', value='20', format='int')
184 op7.addParameter(name='wintitle', value='RTI', format='str')
185 op7.addParameter(name='xmin', value='0', format='float')
186 op7.addParameter(name='xmax', value='24', format='float')
187 op7.addParameter(name='zmin', value='12', format='int')
188 op7.addParameter(name='zmax', value='32', format='int')
189 op7.addParameter(name='showprofile', value='1', format='int')
190 op7.addParameter(name='timerange', value=str(24*60*60), format='int')
191
192 op8 = proc_spectra.addOperation(name='CoherenceMap', optype='external')
193 op8.addParameter(name='id', value='30', format='int')
194 op8.addParameter(name='xmin', value='0.0', format='float')
195 op8.addParameter(name='xmax', value='24.0', format='float')
196
197
198 controller.start()
199
200 if __name__ == '__main__':
201 import time
202 start_time = time.time()
203 main()
204 print("--- %s seconds ---" % (time.time() - start_time))
205
206 '''
207
208
209
210
211
212
213
214
36 215 multiprocess = '''from schainpy.controller import Project, MPProject
37 216 from time import sleep
38 217 desc = "{desc}"
39 218
40 219 ####################
41 220 # PLOTTER RECEIVER #
42 221 ####################
43 222 plotter = Project()
44 223 plotter.setup(id='100', name='receiver', description=desc)
45 224
46 225 receiver_plot = plotter.addProcUnit(name='PlotterReceiver')
47 226 receiver_plot.addParameter(name='throttle', value=20, format='int')
48 227 receiver_plot.addParameter(name='plottypes', value='rti', format='str')
49 228
50 229 rti = receiver_plot.addOperation(name='PlotRTIData', optype='other')
51 230 rti.addParameter(name='zmin', value='-40.0', format='float')
52 231 rti.addParameter(name='zmax', value='100.0', format='float')
53 232 rti.addParameter(name='decimation', value='200', format='int')
54 233 rti.addParameter(name='xmin', value='0.0', format='int')
55 234 rti.addParameter(name='colormap', value='jet', format='str')
56 235
57 236 plotter.start()
58 237
59 238 sleep(2)
60 239
61 240 ################
62 241 # DATA EMITTER #
63 242 ################
64 project = Project()
65 project.setup(id='200', name="{name}", description=desc)
243 controller = Project()
244 controller.setup(id='200', name="{name}", description=desc)
66 245
67 spectra_reader = project.addReadUnit(datatype='SpectraReader',
246 spectra_reader = controller.addReadUnit(datatype='SpectraReader',
68 247 path="{path}",
69 248 startDate={startDate},
70 249 endDate={endDate},
71 250 startTime="{startHour}",
72 251 endTime="{endHour}",
73 252 online=0,
74 253 verbose=1,
75 254 walk=1,
76 255 )
77 256
78 spectra_proc = project.addProcUnit(datatype='Spectra', inputId=spectra_reader.getId())
257 spectra_proc = controller.addProcUnit(datatype='Spectra', inputId=spectra_reader.getId())
79 258
80 parameters_proc = project.addProcUnit(datatype='ParametersProc', inputId=spectra_proc.getId())
259 parameters_proc = controller.addProcUnit(datatype='ParametersProc', inputId=spectra_proc.getId())
81 260 moments = parameters_proc.addOperation(name='SpectralMoments', optype='other')
82 261
83 262 publish = parameters_proc.addOperation(name='PublishData', optype='other')
84 263 publish.addParameter(name='zeromq', value=1, format='int')
85 264 publish.addParameter(name='verbose', value=0, format='bool')
86 265
87 MPProject(project, 16)
266 MPProject(controller, 16)
88 267
89 268
90 269 '''
Show file before | Show file after | Hide comments
@@ -1,1259 +1,1265
1 1 '''
2 2 Updated on January , 2018, for multiprocessing purposes
3 3 Author: Sergio Cortez
4 4 Created on September , 2012
5 5 '''
6 6 from platform import python_version
7 7 import sys
8 8 import ast
9 9 import datetime
10 10 import traceback
11 11 import math
12 12 import time
13 13 import zmq
14 14 from multiprocessing import Process, cpu_count
15 15 from threading import Thread
16 16 from xml.etree.ElementTree import ElementTree, Element, SubElement, tostring
17 17 from xml.dom import minidom
18 18
19 19
20 20 from schainpy.admin import Alarm, SchainWarning
21 21
22 22 ### Temporary imports!!!
23 23 # from schainpy.model import *
24 24 from schainpy.model.io import *
25 25 from schainpy.model.graphics import *
26 26 from schainpy.model.proc.jroproc_base import *
27 27 from schainpy.model.proc.bltrproc_parameters import *
28 28 from schainpy.model.proc.jroproc_spectra import *
29 29 from schainpy.model.proc.jroproc_voltage import *
30 30 from schainpy.model.proc.jroproc_parameters import *
31 31 from schainpy.model.utils.jroutils_publish import *
32 32 from schainpy.utils import log
33 33 ###
34 34
35 35 DTYPES = {
36 36 'Voltage': '.r',
37 37 'Spectra': '.pdata'
38 38 }
39 39
40 40
41 41 def MPProject(project, n=cpu_count()):
42 42 '''
43 43 Project wrapper to run schain in n processes
44 44 '''
45 45
46 46 rconf = project.getReadUnitObj()
47 47 op = rconf.getOperationObj('run')
48 48 dt1 = op.getParameterValue('startDate')
49 49 dt2 = op.getParameterValue('endDate')
50 50 tm1 = op.getParameterValue('startTime')
51 51 tm2 = op.getParameterValue('endTime')
52 52 days = (dt2 - dt1).days
53 53
54 54 for day in range(days + 1):
55 55 skip = 0
56 56 cursor = 0
57 57 processes = []
58 58 dt = dt1 + datetime.timedelta(day)
59 59 dt_str = dt.strftime('%Y/%m/%d')
60 60 reader = JRODataReader()
61 61 paths, files = reader.searchFilesOffLine(path=rconf.path,
62 62 startDate=dt,
63 63 endDate=dt,
64 64 startTime=tm1,
65 65 endTime=tm2,
66 66 ext=DTYPES[rconf.datatype])
67 67 nFiles = len(files)
68 68 if nFiles == 0:
69 69 continue
70 70 skip = int(math.ceil(nFiles / n))
71 71 while nFiles > cursor * skip:
72 72 rconf.update(startDate=dt_str, endDate=dt_str, cursor=cursor,
73 73 skip=skip)
74 74 p = project.clone()
75 75 p.start()
76 76 processes.append(p)
77 77 cursor += 1
78 78
79 79 def beforeExit(exctype, value, trace):
80 80 for process in processes:
81 81 process.terminate()
82 82 process.join()
83 83 print(traceback.print_tb(trace))
84 84
85 85 sys.excepthook = beforeExit
86 86
87 87 for process in processes:
88 88 process.join()
89 89 process.terminate()
90 90
91 91 time.sleep(3)
92 92
93 93 def wait(context):
94 94
95 95 time.sleep(1)
96 96 c = zmq.Context()
97 97 receiver = c.socket(zmq.SUB)
98 98 receiver.connect('ipc:///tmp/schain_{}_pub'.format(self.id))
99 99 receiver.setsockopt(zmq.SUBSCRIBE, self.id.encode())
100 100 log.error('startinggg')
101 101 msg = receiver.recv_multipart()[1]
102 102 #log.error(msg)
103 103 context.terminate()
104 104
105 105 class ParameterConf():
106 106
107 107 id = None
108 108 name = None
109 109 value = None
110 110 format = None
111 111
112 112 __formated_value = None
113 113
114 114 ELEMENTNAME = 'Parameter'
115 115
116 116 def __init__(self):
117 117
118 118 self.format = 'str'
119 119
120 120 def getElementName(self):
121 121
122 122 return self.ELEMENTNAME
123 123
124 124 def getValue(self):
125 125
126 126 value = self.value
127 127 format = self.format
128 128
129 129 if self.__formated_value != None:
130 130
131 131 return self.__formated_value
132 132
133 133 if format == 'obj':
134 134 return value
135 135
136 136 if format == 'str':
137 137 self.__formated_value = str(value)
138 138 return self.__formated_value
139 139
140 140 if value == '':
141 141 raise ValueError('%s: This parameter value is empty' % self.name)
142 142
143 143 if format == 'list':
144 144 strList = value.split(',')
145 145
146 146 self.__formated_value = strList
147 147
148 148 return self.__formated_value
149 149
150 150 if format == 'intlist':
151 151 '''
152 152 Example:
153 153 value = (0,1,2)
154 154 '''
155 155
156 156 new_value = ast.literal_eval(value)
157 157
158 158 if type(new_value) not in (tuple, list):
159 159 new_value = [int(new_value)]
160 160
161 161 self.__formated_value = new_value
162 162
163 163 return self.__formated_value
164 164
165 165 if format == 'floatlist':
166 166 '''
167 167 Example:
168 168 value = (0.5, 1.4, 2.7)
169 169 '''
170 170
171 171 new_value = ast.literal_eval(value)
172 172
173 173 if type(new_value) not in (tuple, list):
174 174 new_value = [float(new_value)]
175 175
176 176 self.__formated_value = new_value
177 177
178 178 return self.__formated_value
179 179
180 180 if format == 'date':
181 181 strList = value.split('/')
182 182 intList = [int(x) for x in strList]
183 183 date = datetime.date(intList[0], intList[1], intList[2])
184 184
185 185 self.__formated_value = date
186 186
187 187 return self.__formated_value
188 188
189 189 if format == 'time':
190 190 strList = value.split(':')
191 191 intList = [int(x) for x in strList]
192 192 time = datetime.time(intList[0], intList[1], intList[2])
193 193
194 194 self.__formated_value = time
195 195
196 196 return self.__formated_value
197 197
198 198 if format == 'pairslist':
199 199 '''
200 200 Example:
201 201 value = (0,1),(1,2)
202 202 '''
203 203
204 204 new_value = ast.literal_eval(value)
205 205
206 206 if type(new_value) not in (tuple, list):
207 207 raise ValueError('%s has to be a tuple or list of pairs' % value)
208 208
209 209 if type(new_value[0]) not in (tuple, list):
210 210 if len(new_value) != 2:
211 211 raise ValueError('%s has to be a tuple or list of pairs' % value)
212 212 new_value = [new_value]
213 213
214 214 for thisPair in new_value:
215 215 if len(thisPair) != 2:
216 216 raise ValueError('%s has to be a tuple or list of pairs' % value)
217 217
218 218 self.__formated_value = new_value
219 219
220 220 return self.__formated_value
221 221
222 222 if format == 'multilist':
223 223 '''
224 224 Example:
225 225 value = (0,1,2),(3,4,5)
226 226 '''
227 227 multiList = ast.literal_eval(value)
228 228
229 229 if type(multiList[0]) == int:
230 230 multiList = ast.literal_eval('(' + value + ')')
231 231
232 232 self.__formated_value = multiList
233 233
234 234 return self.__formated_value
235 235
236 236 if format == 'bool':
237 237 value = int(value)
238 238
239 239 if format == 'int':
240 240 value = float(value)
241 241
242 242 format_func = eval(format)
243 243
244 244 self.__formated_value = format_func(value)
245 245
246 246 return self.__formated_value
247 247
248 248 def updateId(self, new_id):
249 249
250 250 self.id = str(new_id)
251 251
252 252 def setup(self, id, name, value, format='str'):
253 253 self.id = str(id)
254 254 self.name = name
255 255 if format == 'obj':
256 256 self.value = value
257 257 else:
258 258 self.value = str(value)
259 259 self.format = str.lower(format)
260 260
261 261 self.getValue()
262 262
263 263 return 1
264 264
265 265 def update(self, name, value, format='str'):
266 266
267 267 self.name = name
268 268 self.value = str(value)
269 269 self.format = format
270 270
271 271 def makeXml(self, opElement):
272 272 if self.name not in ('queue',):
273 273 parmElement = SubElement(opElement, self.ELEMENTNAME)
274 274 parmElement.set('id', str(self.id))
275 275 parmElement.set('name', self.name)
276 276 parmElement.set('value', self.value)
277 277 parmElement.set('format', self.format)
278
278
279 279 def readXml(self, parmElement):
280 280
281 281 self.id = parmElement.get('id')
282 282 self.name = parmElement.get('name')
283 283 self.value = parmElement.get('value')
284 284 self.format = str.lower(parmElement.get('format'))
285 285
286 286 # Compatible with old signal chain version
287 287 if self.format == 'int' and self.name == 'idfigure':
288 288 self.name = 'id'
289 289
290 290 def printattr(self):
291 291
292 print('Parameter[%s]: name = %s, value = %s, format = %s' % (self.id, self.name, self.value, self.format))
292 print('Parameter[%s]: name = %s, value = %s, format = %s, project_id = %s' % (self.id, self.name, self.value, self.format, self.project_id))
293 293
294 294 class OperationConf():
295 295
296 296 ELEMENTNAME = 'Operation'
297 297
298 298 def __init__(self):
299 299
300 300 self.id = '0'
301 301 self.name = None
302 302 self.priority = None
303 303 self.topic = None
304 304
305 305 def __getNewId(self):
306 306
307 307 return int(self.id) * 10 + len(self.parmConfObjList) + 1
308 308
309 309 def getId(self):
310 310 return self.id
311 311
312 312 def updateId(self, new_id):
313 313
314 314 self.id = str(new_id)
315 315
316 316 n = 1
317 317 for parmObj in self.parmConfObjList:
318 318
319 319 idParm = str(int(new_id) * 10 + n)
320 320 parmObj.updateId(idParm)
321 321
322 322 n += 1
323 323
324 324 def getElementName(self):
325 325
326 326 return self.ELEMENTNAME
327 327
328 328 def getParameterObjList(self):
329 329
330 330 return self.parmConfObjList
331 331
332 332 def getParameterObj(self, parameterName):
333 333
334 334 for parmConfObj in self.parmConfObjList:
335 335
336 336 if parmConfObj.name != parameterName:
337 337 continue
338 338
339 339 return parmConfObj
340 340
341 341 return None
342 342
343 343 def getParameterObjfromValue(self, parameterValue):
344 344
345 345 for parmConfObj in self.parmConfObjList:
346 346
347 347 if parmConfObj.getValue() != parameterValue:
348 348 continue
349 349
350 350 return parmConfObj.getValue()
351 351
352 352 return None
353 353
354 354 def getParameterValue(self, parameterName):
355 355
356 356 parameterObj = self.getParameterObj(parameterName)
357 357
358 358 # if not parameterObj:
359 359 # return None
360 360
361 361 value = parameterObj.getValue()
362 362
363 363 return value
364 364
365 365 def getKwargs(self):
366 366
367 367 kwargs = {}
368 368
369 369 for parmConfObj in self.parmConfObjList:
370 370 if self.name == 'run' and parmConfObj.name == 'datatype':
371 371 continue
372 372
373 373 kwargs[parmConfObj.name] = parmConfObj.getValue()
374 374
375 375 return kwargs
376 376
377 377 def setup(self, id, name, priority, type, project_id):
378 378
379 379 self.id = str(id)
380 380 self.project_id = project_id
381 381 self.name = name
382 382 self.type = type
383 383 self.priority = priority
384 384 self.parmConfObjList = []
385 385
386 386 def removeParameters(self):
387 387
388 388 for obj in self.parmConfObjList:
389 389 del obj
390 390
391 391 self.parmConfObjList = []
392 392
393 393 def addParameter(self, name, value, format='str'):
394 394
395 395 if value is None:
396 396 return None
397 397 id = self.__getNewId()
398 398
399 399 parmConfObj = ParameterConf()
400 400 if not parmConfObj.setup(id, name, value, format):
401 401 return None
402 402
403 403 self.parmConfObjList.append(parmConfObj)
404 404
405 405 return parmConfObj
406 406
407 407 def changeParameter(self, name, value, format='str'):
408 408
409 409 parmConfObj = self.getParameterObj(name)
410 410 parmConfObj.update(name, value, format)
411 411
412 412 return parmConfObj
413 413
414 414 def makeXml(self, procUnitElement):
415 415
416 416 opElement = SubElement(procUnitElement, self.ELEMENTNAME)
417 417 opElement.set('id', str(self.id))
418 418 opElement.set('name', self.name)
419 419 opElement.set('type', self.type)
420 420 opElement.set('priority', str(self.priority))
421 421
422 422 for parmConfObj in self.parmConfObjList:
423 423 parmConfObj.makeXml(opElement)
424 424
425 def readXml(self, opElement):
425 def readXml(self, opElement, project_id):
426 426
427 427 self.id = opElement.get('id')
428 428 self.name = opElement.get('name')
429 429 self.type = opElement.get('type')
430 430 self.priority = opElement.get('priority')
431 self.project_id = str(project_id) #yong
431 432
432 433 # Compatible with old signal chain version
433 434 # Use of 'run' method instead 'init'
434 435 if self.type == 'self' and self.name == 'init':
435 436 self.name = 'run'
436 437
437 438 self.parmConfObjList = []
438 439
439 440 parmElementList = opElement.iter(ParameterConf().getElementName())
440 441
441 442 for parmElement in parmElementList:
442 443 parmConfObj = ParameterConf()
443 444 parmConfObj.readXml(parmElement)
444 445
445 446 # Compatible with old signal chain version
446 447 # If an 'plot' OPERATION is found, changes name operation by the value of its type PARAMETER
447 448 if self.type != 'self' and self.name == 'Plot':
448 449 if parmConfObj.format == 'str' and parmConfObj.name == 'type':
449 450 self.name = parmConfObj.value
450 451 continue
451 452
452 453 self.parmConfObjList.append(parmConfObj)
453 454
454 455 def printattr(self):
455 456
456 print('%s[%s]: name = %s, type = %s, priority = %s' % (self.ELEMENTNAME,
457 print('%s[%s]: name = %s, type = %s, priority = %s, project_id = %s' % (self.ELEMENTNAME,
457 458 self.id,
458 459 self.name,
459 460 self.type,
460 self.priority))
461 self.priority,
462 self.project_id))
461 463
462 464 for parmConfObj in self.parmConfObjList:
463 465 parmConfObj.printattr()
464 466
465 467 def createObject(self):
466 468
467 469 className = eval(self.name)
468
470
469 471 if self.type == 'other':
470 472 opObj = className()
471 473 elif self.type == 'external':
472 474 kwargs = self.getKwargs()
473 475 opObj = className(self.id, self.project_id, **kwargs)
474 476 opObj.start()
475 477
476 478 return opObj
477 479
478 480 class ProcUnitConf():
479 481
480 482 ELEMENTNAME = 'ProcUnit'
481 483
482 484 def __init__(self):
483 485
484 486 self.id = None
485 487 self.datatype = None
486 488 self.name = None
487 489 self.inputId = None
488 490 self.opConfObjList = []
489 491 self.procUnitObj = None
490 492 self.opObjDict = {}
491 493
492 494 def __getPriority(self):
493 495
494 496 return len(self.opConfObjList) + 1
495 497
496 498 def __getNewId(self):
497 499
498 500 return int(self.id) * 10 + len(self.opConfObjList) + 1
499 501
500 502 def getElementName(self):
501 503
502 504 return self.ELEMENTNAME
503 505
504 506 def getId(self):
505 507
506 508 return self.id
507 509
508 510 def updateId(self, new_id):
509 511 '''
510 512 new_id = int(parentId) * 10 + (int(self.id) % 10)
511 513 new_inputId = int(parentId) * 10 + (int(self.inputId) % 10)
512 514
513 515 # If this proc unit has not inputs
514 516 #if self.inputId == '0':
515 517 #new_inputId = 0
516 518
517 519 n = 1
518 520 for opConfObj in self.opConfObjList:
519 521
520 522 idOp = str(int(new_id) * 10 + n)
521 523 opConfObj.updateId(idOp)
522 524
523 525 n += 1
524 526
525 527 self.parentId = str(parentId)
526 528 self.id = str(new_id)
527 529 #self.inputId = str(new_inputId)
528 530 '''
529 531 n = 1
530 532
531 533 def getInputId(self):
532 534
533 535 return self.inputId
534 536
535 537 def getOperationObjList(self):
536 538
537 539 return self.opConfObjList
538 540
539 541 def getOperationObj(self, name=None):
540 542
541 543 for opConfObj in self.opConfObjList:
542 544
543 545 if opConfObj.name != name:
544 546 continue
545 547
546 548 return opConfObj
547 549
548 550 return None
549 551
550 552 def getOpObjfromParamValue(self, value=None):
551 553
552 554 for opConfObj in self.opConfObjList:
553 555 if opConfObj.getParameterObjfromValue(parameterValue=value) != value:
554 556 continue
555 557 return opConfObj
556 558 return None
557 559
558 560 def getProcUnitObj(self):
559 561
560 562 return self.procUnitObj
561 563
562 564 def setup(self, project_id, id, name, datatype, inputId):
563 565 '''
564 566 id sera el topico a publicar
565 567 inputId sera el topico a subscribirse
566 568 '''
567 569
568 570 # Compatible with old signal chain version
569 571 if datatype == None and name == None:
570 572 raise ValueError('datatype or name should be defined')
571 573
572 574 #Definir una condicion para inputId cuando sea 0
573 575
574 576 if name == None:
575 577 if 'Proc' in datatype:
576 578 name = datatype
577 579 else:
578 580 name = '%sProc' % (datatype)
579 581
580 582 if datatype == None:
581 583 datatype = name.replace('Proc', '')
582 584
583 585 self.id = str(id)
584 586 self.project_id = project_id
585 587 self.name = name
586 588 self.datatype = datatype
587 589 self.inputId = inputId
588 590 self.opConfObjList = []
589 591
590 592 self.addOperation(name='run', optype='self')
591 593
592 594 def removeOperations(self):
593 595
594 596 for obj in self.opConfObjList:
595 597 del obj
596 598
597 599 self.opConfObjList = []
598 600 self.addOperation(name='run')
599 601
600 602 def addParameter(self, **kwargs):
601 603 '''
602 604 Add parameters to 'run' operation
603 605 '''
604 606 opObj = self.opConfObjList[0]
605 607
606 608 opObj.addParameter(**kwargs)
607 609
608 610 return opObj
609 611
610 612 def addOperation(self, name, optype='self'):
611 613 '''
612 614 Actualizacion - > proceso comunicacion
613 615 En el caso de optype='self', elminar. DEfinir comuncacion IPC -> Topic
614 616 definir el tipoc de socket o comunicacion ipc++
615 617
616 618 '''
617 619
618 620 id = self.__getNewId()
619 621 priority = self.__getPriority() # Sin mucho sentido, pero puede usarse
620 622 opConfObj = OperationConf()
621 623 opConfObj.setup(id, name=name, priority=priority, type=optype, project_id=self.project_id)
622 624 self.opConfObjList.append(opConfObj)
623 625
624 626 return opConfObj
625 627
626 628 def makeXml(self, projectElement):
627 629
628 630 procUnitElement = SubElement(projectElement, self.ELEMENTNAME)
629 631 procUnitElement.set('id', str(self.id))
630 632 procUnitElement.set('name', self.name)
631 633 procUnitElement.set('datatype', self.datatype)
632 634 procUnitElement.set('inputId', str(self.inputId))
633 635
634 636 for opConfObj in self.opConfObjList:
635 637 opConfObj.makeXml(procUnitElement)
636 638
637 def readXml(self, upElement):
639 def readXml(self, upElement, project_id):
638 640
639 641 self.id = upElement.get('id')
640 642 self.name = upElement.get('name')
641 643 self.datatype = upElement.get('datatype')
642 644 self.inputId = upElement.get('inputId')
645 self.project_id = str(project_id)
643 646
644 647 if self.ELEMENTNAME == 'ReadUnit':
645 648 self.datatype = self.datatype.replace('Reader', '')
646 649
647 650 if self.ELEMENTNAME == 'ProcUnit':
648 651 self.datatype = self.datatype.replace('Proc', '')
649 652
650 653 if self.inputId == 'None':
651 654 self.inputId = '0'
652 655
653 656 self.opConfObjList = []
654 657
655 658 opElementList = upElement.iter(OperationConf().getElementName())
656 659
657 660 for opElement in opElementList:
658 661 opConfObj = OperationConf()
659 opConfObj.readXml(opElement)
662 opConfObj.readXml(opElement, project_id)
660 663 self.opConfObjList.append(opConfObj)
661 664
662 665 def printattr(self):
663 666
664 print('%s[%s]: name = %s, datatype = %s, inputId = %s' % (self.ELEMENTNAME,
667 print('%s[%s]: name = %s, datatype = %s, inputId = %s, project_id = %s' % (self.ELEMENTNAME,
665 668 self.id,
666 669 self.name,
667 670 self.datatype,
668 self.inputId))
671 self.inputId,
672 self.project_id))
669 673
670 674 for opConfObj in self.opConfObjList:
671 675 opConfObj.printattr()
672 676
673 677 def getKwargs(self):
674 678
675 679 opObj = self.opConfObjList[0]
676 680 kwargs = opObj.getKwargs()
677 681
678 682 return kwargs
679 683
680 684 def createObjects(self):
681 685 '''
682 686 Instancia de unidades de procesamiento.
683 687 '''
684 688 className = eval(self.name)
685 689 kwargs = self.getKwargs()
686 690 procUnitObj = className(self.id, self.inputId, self.project_id, **kwargs) # necesitan saber su id y su entrada por fines de ipc
687 691 log.success('creating process...', self.name)
688 692
689 693 for opConfObj in self.opConfObjList:
690 694
691 695 if opConfObj.type == 'self' and opConfObj.name == 'run':
692 696 continue
693 697 elif opConfObj.type == 'self':
694 698 opObj = getattr(procUnitObj, opConfObj.name)
695 699 else:
696 700 opObj = opConfObj.createObject()
697 701
698 702 log.success('creating operation: {}, type:{}'.format(
699 703 opConfObj.name,
700 704 opConfObj.type), self.name)
701 705
702 706 procUnitObj.addOperation(opConfObj, opObj)
703 707
704 708 procUnitObj.start()
705 709 self.procUnitObj = procUnitObj
706 710
707 711 def close(self):
708 712
709 713 for opConfObj in self.opConfObjList:
710 714 if opConfObj.type == 'self':
711 715 continue
712 716
713 717 opObj = self.procUnitObj.getOperationObj(opConfObj.id)
714 718 opObj.close()
715 719
716 720 self.procUnitObj.close()
717 721
718 722 return
719 723
720 724
721 725 class ReadUnitConf(ProcUnitConf):
722 726
723 727 ELEMENTNAME = 'ReadUnit'
724 728
725 729 def __init__(self):
726 730
727 731 self.id = None
728 732 self.datatype = None
729 733 self.name = None
730 734 self.inputId = None
731 735 self.opConfObjList = []
732 736
733 737 def getElementName(self):
734 738
735 739 return self.ELEMENTNAME
736 740
737 741 def setup(self, project_id, id, name, datatype, path='', startDate='', endDate='',
738 742 startTime='', endTime='', server=None, **kwargs):
739 743
740 744
741 745 '''
742 746 *****el id del proceso sera el Topico
743 747
744 748 Adicion de {topic}, si no esta presente -> error
745 749 kwargs deben ser trasmitidos en la instanciacion
746 750
747 751 '''
748 752
749 753 # Compatible with old signal chain version
750 754 if datatype == None and name == None:
751 755 raise ValueError('datatype or name should be defined')
752 756 if name == None:
753 757 if 'Reader' in datatype:
754 758 name = datatype
755 759 datatype = name.replace('Reader','')
756 760 else:
757 761 name = '{}Reader'.format(datatype)
758 762 if datatype == None:
759 763 if 'Reader' in name:
760 764 datatype = name.replace('Reader','')
761 765 else:
762 766 datatype = name
763 767 name = '{}Reader'.format(name)
764 768
765 769 self.id = id
766 770 self.project_id = project_id
767 771 self.name = name
768 772 self.datatype = datatype
769 773 if path != '':
770 774 self.path = os.path.abspath(path)
771 775 self.startDate = startDate
772 776 self.endDate = endDate
773 777 self.startTime = startTime
774 778 self.endTime = endTime
775 779 self.server = server
776 780 self.addRunOperation(**kwargs)
777 781
778 782 def update(self, **kwargs):
779 783
780 784 if 'datatype' in kwargs:
781 785 datatype = kwargs.pop('datatype')
782 786 if 'Reader' in datatype:
783 787 self.name = datatype
784 788 else:
785 789 self.name = '%sReader' % (datatype)
786 790 self.datatype = self.name.replace('Reader', '')
787 791
788 792 attrs = ('path', 'startDate', 'endDate',
789 793 'startTime', 'endTime')
790 794
791 795 for attr in attrs:
792 796 if attr in kwargs:
793 797 setattr(self, attr, kwargs.pop(attr))
794 798
795 799 self.updateRunOperation(**kwargs)
796 800
797 801 def removeOperations(self):
798 802
799 803 for obj in self.opConfObjList:
800 804 del obj
801 805
802 806 self.opConfObjList = []
803 807
804 808 def addRunOperation(self, **kwargs):
805 809
806 810 opObj = self.addOperation(name='run', optype='self')
807 811
808 812 if self.server is None:
809 813 opObj.addParameter(
810 814 name='datatype', value=self.datatype, format='str')
811 815 opObj.addParameter(name='path', value=self.path, format='str')
812 816 opObj.addParameter(
813 817 name='startDate', value=self.startDate, format='date')
814 818 opObj.addParameter(
815 819 name='endDate', value=self.endDate, format='date')
816 820 opObj.addParameter(
817 821 name='startTime', value=self.startTime, format='time')
818 822 opObj.addParameter(
819 823 name='endTime', value=self.endTime, format='time')
820 824
821 825 for key, value in list(kwargs.items()):
822 826 opObj.addParameter(name=key, value=value,
823 827 format=type(value).__name__)
824 828 else:
825 829 opObj.addParameter(name='server', value=self.server, format='str')
826 830
827 831 return opObj
828 832
829 833 def updateRunOperation(self, **kwargs):
830 834
831 835 opObj = self.getOperationObj(name='run')
832 836 opObj.removeParameters()
833 837
834 838 opObj.addParameter(name='datatype', value=self.datatype, format='str')
835 839 opObj.addParameter(name='path', value=self.path, format='str')
836 840 opObj.addParameter(
837 841 name='startDate', value=self.startDate, format='date')
838 842 opObj.addParameter(name='endDate', value=self.endDate, format='date')
839 843 opObj.addParameter(
840 844 name='startTime', value=self.startTime, format='time')
841 845 opObj.addParameter(name='endTime', value=self.endTime, format='time')
842 846
843 847 for key, value in list(kwargs.items()):
844 848 opObj.addParameter(name=key, value=value,
845 849 format=type(value).__name__)
846 850
847 851 return opObj
848 852
849 def readXml(self, upElement):
853 def readXml(self, upElement, project_id):
850 854
851 855 self.id = upElement.get('id')
852 856 self.name = upElement.get('name')
853 857 self.datatype = upElement.get('datatype')
858 self.project_id = str(project_id) #yong
854 859
855 860 if self.ELEMENTNAME == 'ReadUnit':
856 861 self.datatype = self.datatype.replace('Reader', '')
857 862
858 863 self.opConfObjList = []
859 864
860 865 opElementList = upElement.iter(OperationConf().getElementName())
861 866
862 867 for opElement in opElementList:
863 868 opConfObj = OperationConf()
864 opConfObj.readXml(opElement)
869 opConfObj.readXml(opElement, project_id)
865 870 self.opConfObjList.append(opConfObj)
866 871
867 872 if opConfObj.name == 'run':
868 873 self.path = opConfObj.getParameterValue('path')
869 874 self.startDate = opConfObj.getParameterValue('startDate')
870 875 self.endDate = opConfObj.getParameterValue('endDate')
871 876 self.startTime = opConfObj.getParameterValue('startTime')
872 877 self.endTime = opConfObj.getParameterValue('endTime')
873 878
874 879
875 880 class Project(Process):
876 881
877 882 ELEMENTNAME = 'Project'
878 883
879 884 def __init__(self):
880 885
881 886 Process.__init__(self)
882 887 self.id = None
883 888 self.filename = None
884 889 self.description = None
885 890 self.email = None
886 891 self.alarm = None
887 892 self.procUnitConfObjDict = {}
888 893
889 894 def __getNewId(self):
890 895
891 896 idList = list(self.procUnitConfObjDict.keys())
892 897 id = int(self.id) * 10
893 898
894 899 while True:
895 900 id += 1
896 901
897 902 if str(id) in idList:
898 903 continue
899 904
900 905 break
901 906
902 907 return str(id)
903 908
904 909 def getElementName(self):
905 910
906 911 return self.ELEMENTNAME
907 912
908 913 def getId(self):
909 914
910 915 return self.id
911 916
912 917 def updateId(self, new_id):
913 918
914 919 self.id = str(new_id)
915 920
916 921 keyList = list(self.procUnitConfObjDict.keys())
917 922 keyList.sort()
918 923
919 924 n = 1
920 925 newProcUnitConfObjDict = {}
921 926
922 927 for procKey in keyList:
923 928
924 929 procUnitConfObj = self.procUnitConfObjDict[procKey]
925 930 idProcUnit = str(int(self.id) * 10 + n)
926 931 procUnitConfObj.updateId(idProcUnit)
927 932 newProcUnitConfObjDict[idProcUnit] = procUnitConfObj
928 933 n += 1
929 934
930 935 self.procUnitConfObjDict = newProcUnitConfObjDict
931 936
932 937 def setup(self, id=1, name='', description='', email=None, alarm=[]):
933 938
934 939 print(' ')
935 940 print('*' * 60)
936 941 print('* Starting SIGNAL CHAIN PROCESSING (Multiprocessing) v%s *' % schainpy.__version__)
937 942 print('*' * 60)
938 943 print("* Python " + python_version() + " *")
939 944 print('*' * 19)
940 945 print(' ')
941 946 self.id = str(id)
942 947 self.description = description
943 948 self.email = email
944 949 self.alarm = alarm
945 950
946 951 def update(self, **kwargs):
947 952
948 953 for key, value in list(kwargs.items()):
949 954 setattr(self, key, value)
950 955
951 956 def clone(self):
952 957
953 958 p = Project()
954 959 p.procUnitConfObjDict = self.procUnitConfObjDict
955 960 return p
956 961
957 962 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
958 963
959 964 '''
960 965 Actualizacion:
961 966 Se agrego un nuevo argumento: topic -relativo a la forma de comunicar los procesos simultaneos
962 967
963 968 * El id del proceso sera el topico al que se deben subscribir los procUnits para recibir la informacion(data)
964 969
965 970 '''
966 971
967 972 if id is None:
968 973 idReadUnit = self.__getNewId()
969 974 else:
970 975 idReadUnit = str(id)
971 976
972 977 readUnitConfObj = ReadUnitConf()
973 978 readUnitConfObj.setup(self.id, idReadUnit, name, datatype, **kwargs)
974 979 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
975 980
976 981 return readUnitConfObj
977 982
978 983 def addProcUnit(self, inputId='0', datatype=None, name=None):
979 984
980 985 '''
981 986 Actualizacion:
982 987 Se agrego dos nuevos argumentos: topic_read (lee data de otro procUnit) y topic_write(escribe o envia data a otro procUnit)
983 988 Deberia reemplazar a "inputId"
984 989
985 990 ** A fin de mantener el inputID, este sera la representaacion del topicoal que deben subscribirse. El ID propio de la intancia
986 991 (proceso) sera el topico de la publicacion, todo sera asignado de manera dinamica.
987 992
988 993 '''
989 994
990 995 idProcUnit = self.__getNewId() #Topico para subscripcion
991 996 procUnitConfObj = ProcUnitConf()
992 997 procUnitConfObj.setup(self.id, idProcUnit, name, datatype, inputId) #topic_read, topic_write,
993 998 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
994 999
995 1000 return procUnitConfObj
996 1001
997 1002 def removeProcUnit(self, id):
998 1003
999 1004 if id in list(self.procUnitConfObjDict.keys()):
1000 1005 self.procUnitConfObjDict.pop(id)
1001 1006
1002 1007 def getReadUnitId(self):
1003 1008
1004 1009 readUnitConfObj = self.getReadUnitObj()
1005 1010
1006 1011 return readUnitConfObj.id
1007 1012
1008 1013 def getReadUnitObj(self):
1009 1014
1010 1015 for obj in list(self.procUnitConfObjDict.values()):
1011 1016 if obj.getElementName() == 'ReadUnit':
1012 1017 return obj
1013 1018
1014 1019 return None
1015 1020
1016 1021 def getProcUnitObj(self, id=None, name=None):
1017 1022
1018 1023 if id != None:
1019 1024 return self.procUnitConfObjDict[id]
1020 1025
1021 1026 if name != None:
1022 1027 return self.getProcUnitObjByName(name)
1023 1028
1024 1029 return None
1025 1030
1026 1031 def getProcUnitObjByName(self, name):
1027 1032
1028 1033 for obj in list(self.procUnitConfObjDict.values()):
1029 1034 if obj.name == name:
1030 1035 return obj
1031 1036
1032 1037 return None
1033 1038
1034 1039 def procUnitItems(self):
1035 1040
1036 1041 return list(self.procUnitConfObjDict.items())
1037 1042
1038 1043 def makeXml(self):
1039 1044
1040 1045 projectElement = Element('Project')
1041 1046 projectElement.set('id', str(self.id))
1042 1047 projectElement.set('name', self.name)
1043 1048 projectElement.set('description', self.description)
1044 1049
1045 1050 for procUnitConfObj in list(self.procUnitConfObjDict.values()):
1046 1051 procUnitConfObj.makeXml(projectElement)
1047 1052
1048 1053 self.projectElement = projectElement
1049 1054
1050 1055 def writeXml(self, filename=None):
1051 1056
1052 1057 if filename == None:
1053 1058 if self.filename:
1054 1059 filename = self.filename
1055 1060 else:
1056 1061 filename = 'schain.xml'
1057 1062
1058 1063 if not filename:
1059 1064 print('filename has not been defined. Use setFilename(filename) for do it.')
1060 1065 return 0
1061 1066
1062 1067 abs_file = os.path.abspath(filename)
1063 1068
1064 1069 if not os.access(os.path.dirname(abs_file), os.W_OK):
1065 1070 print('No write permission on %s' % os.path.dirname(abs_file))
1066 1071 return 0
1067 1072
1068 1073 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
1069 1074 print('File %s already exists and it could not be overwriten' % abs_file)
1070 1075 return 0
1071 1076
1072 1077 self.makeXml()
1073 1078
1074 1079 ElementTree(self.projectElement).write(abs_file, method='xml')
1075 1080
1076 1081 self.filename = abs_file
1077 1082
1078 1083 return 1
1079 1084
1080 1085 def readXml(self, filename=None):
1081 1086
1082 1087 if not filename:
1083 1088 print('filename is not defined')
1084 1089 return 0
1085 1090
1086 1091 abs_file = os.path.abspath(filename)
1087 1092
1088 1093 if not os.path.isfile(abs_file):
1089 1094 print('%s file does not exist' % abs_file)
1090 1095 return 0
1091 1096
1092 1097 self.projectElement = None
1093 1098 self.procUnitConfObjDict = {}
1094 1099
1095 1100 try:
1096 1101 self.projectElement = ElementTree().parse(abs_file)
1097 1102 except:
1098 1103 print('Error reading %s, verify file format' % filename)
1099 1104 return 0
1100 1105
1101 1106 self.project = self.projectElement.tag
1102 1107
1103 1108 self.id = self.projectElement.get('id')
1104 1109 self.name = self.projectElement.get('name')
1105 1110 self.description = self.projectElement.get('description')
1106 1111
1107 1112 readUnitElementList = self.projectElement.iter(
1108 1113 ReadUnitConf().getElementName())
1109 1114
1110 1115 for readUnitElement in readUnitElementList:
1111 1116 readUnitConfObj = ReadUnitConf()
1112 readUnitConfObj.readXml(readUnitElement)
1117 readUnitConfObj.readXml(readUnitElement, self.id)
1113 1118 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
1114 1119
1115 1120 procUnitElementList = self.projectElement.iter(
1116 1121 ProcUnitConf().getElementName())
1117 1122
1118 1123 for procUnitElement in procUnitElementList:
1119 1124 procUnitConfObj = ProcUnitConf()
1120 procUnitConfObj.readXml(procUnitElement)
1125 procUnitConfObj.readXml(procUnitElement, self.id)
1121 1126 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
1122 1127
1123 1128 self.filename = abs_file
1124 1129
1125 1130 return 1
1126 1131
1127 1132 def __str__(self):
1128 1133
1129 print('Project[%s]: name = %s, description = %s' % (self.id,
1134 print('Project[%s]: name = %s, description = %s, project_id = %s' % (self.id,
1130 1135 self.name,
1131 self.description))
1136 self.description,
1137 self.project_id))
1132 1138
1133 1139 for procUnitConfObj in self.procUnitConfObjDict.values():
1134 1140 print(procUnitConfObj)
1135 1141
1136 1142 def createObjects(self):
1137 1143
1138 1144 for procUnitConfObj in self.procUnitConfObjDict.values():
1139 1145 procUnitConfObj.createObjects()
1140 1146
1141 1147 def __handleError(self, procUnitConfObj, modes=None, stdout=True):
1142 1148
1143 1149 import socket
1144 1150
1145 1151 if modes is None:
1146 1152 modes = self.alarm
1147 1153
1148 1154 if not self.alarm:
1149 1155 modes = []
1150 1156
1151 1157 err = traceback.format_exception(sys.exc_info()[0],
1152 1158 sys.exc_info()[1],
1153 1159 sys.exc_info()[2])
1154 1160
1155 1161 log.error('{}'.format(err[-1]), procUnitConfObj.name)
1156 1162
1157 1163 message = ''.join(err)
1158 1164
1159 1165 if stdout:
1160 1166 sys.stderr.write(message)
1161 1167
1162 1168 subject = 'SChain v%s: Error running %s\n' % (
1163 1169 schainpy.__version__, procUnitConfObj.name)
1164 1170
1165 1171 subtitle = '%s: %s\n' % (
1166 1172 procUnitConfObj.getElementName(), procUnitConfObj.name)
1167 1173 subtitle += 'Hostname: %s\n' % socket.gethostbyname(
1168 1174 socket.gethostname())
1169 1175 subtitle += 'Working directory: %s\n' % os.path.abspath('./')
1170 1176 subtitle += 'Configuration file: %s\n' % self.filename
1171 1177 subtitle += 'Time: %s\n' % str(datetime.datetime.now())
1172 1178
1173 1179 readUnitConfObj = self.getReadUnitObj()
1174 1180 if readUnitConfObj:
1175 1181 subtitle += '\nInput parameters:\n'
1176 1182 subtitle += '[Data path = %s]\n' % readUnitConfObj.path
1177 1183 subtitle += '[Data type = %s]\n' % readUnitConfObj.datatype
1178 1184 subtitle += '[Start date = %s]\n' % readUnitConfObj.startDate
1179 1185 subtitle += '[End date = %s]\n' % readUnitConfObj.endDate
1180 1186 subtitle += '[Start time = %s]\n' % readUnitConfObj.startTime
1181 1187 subtitle += '[End time = %s]\n' % readUnitConfObj.endTime
1182 1188
1183 1189 a = Alarm(
1184 1190 modes=modes,
1185 1191 email=self.email,
1186 1192 message=message,
1187 1193 subject=subject,
1188 1194 subtitle=subtitle,
1189 1195 filename=self.filename
1190 1196 )
1191 1197
1192 1198 return a
1193 1199
1194 1200 def isPaused(self):
1195 1201 return 0
1196 1202
1197 1203 def isStopped(self):
1198 1204 return 0
1199 1205
1200 1206 def runController(self):
1201 1207 '''
1202 1208 returns 0 when this process has been stopped, 1 otherwise
1203 1209 '''
1204 1210
1205 1211 if self.isPaused():
1206 1212 print('Process suspended')
1207 1213
1208 1214 while True:
1209 1215 time.sleep(0.1)
1210 1216
1211 1217 if not self.isPaused():
1212 1218 break
1213 1219
1214 1220 if self.isStopped():
1215 1221 break
1216 1222
1217 1223 print('Process reinitialized')
1218 1224
1219 1225 if self.isStopped():
1220 1226 print('Process stopped')
1221 1227 return 0
1222 1228
1223 1229 return 1
1224 1230
1225 1231 def setFilename(self, filename):
1226 1232
1227 1233 self.filename = filename
1228 1234
1229 1235 def setProxyCom(self):
1230 1236
1231 1237 if not os.path.exists('/tmp/schain'):
1232 1238 os.mkdir('/tmp/schain')
1233 1239
1234 1240 self.ctx = zmq.Context()
1235 1241 xpub = self.ctx.socket(zmq.XPUB)
1236 1242 xpub.bind('ipc:///tmp/schain/{}_pub'.format(self.id))
1237 1243 xsub = self.ctx.socket(zmq.XSUB)
1238 1244 xsub.bind('ipc:///tmp/schain/{}_sub'.format(self.id))
1239 1245
1240 1246 try:
1241 1247 zmq.proxy(xpub, xsub)
1242 1248 except zmq.ContextTerminated:
1243 1249 xpub.close()
1244 1250 xsub.close()
1245 1251
1246 1252 def run(self):
1247 1253
1248 1254 log.success('Starting {}: {}'.format(self.name, self.id), tag='')
1249 1255 self.start_time = time.time()
1250 1256 self.createObjects()
1251 1257 # t = Thread(target=wait, args=(self.ctx, ))
1252 1258 # t.start()
1253 1259 self.setProxyCom()
1254 1260
1255 1261 # Iniciar todos los procesos .start(), monitoreo de procesos. ELiminar lo de abajo
1256 1262
1257 1263 log.success('{} finished (time: {}s)'.format(
1258 1264 self.name,
1259 1265 time.time()-self.start_time))
Show file before | Show file after | Hide comments
@@ -1,367 +1,369
1 1 '''
2 2 Created on Nov 9, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7
8 8 import os
9 9 import sys
10 10 import time
11 11 import glob
12 12 import datetime
13 13
14 14 import numpy
15 15
16 16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
17 17 from schainpy.model.data.jrodata import Parameters
18 18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 19 from schainpy.utils import log
20 20
21 21 FILE_HEADER_STRUCTURE = numpy.dtype([
22 22 ('FMN', '<u4'),
23 23 ('nrec', '<u4'),
24 24 ('fr_offset', '<u4'),
25 25 ('id', '<u4'),
26 26 ('site', 'u1', (32,))
27 27 ])
28 28
29 29 REC_HEADER_STRUCTURE = numpy.dtype([
30 30 ('rmn', '<u4'),
31 31 ('rcounter', '<u4'),
32 32 ('nr_offset', '<u4'),
33 33 ('tr_offset', '<u4'),
34 34 ('time', '<u4'),
35 35 ('time_msec', '<u4'),
36 36 ('tag', 'u1', (32,)),
37 37 ('comments', 'u1', (32,)),
38 38 ('lat', '<f4'),
39 39 ('lon', '<f4'),
40 40 ('gps_status', '<u4'),
41 41 ('freq', '<u4'),
42 42 ('freq0', '<u4'),
43 43 ('nchan', '<u4'),
44 44 ('delta_r', '<u4'),
45 45 ('nranges', '<u4'),
46 46 ('r0', '<u4'),
47 47 ('prf', '<u4'),
48 48 ('ncoh', '<u4'),
49 49 ('npoints', '<u4'),
50 50 ('polarization', '<i4'),
51 51 ('rx_filter', '<u4'),
52 52 ('nmodes', '<u4'),
53 53 ('dmode_index', '<u4'),
54 54 ('dmode_rngcorr', '<u4'),
55 55 ('nrxs', '<u4'),
56 56 ('acf_length', '<u4'),
57 57 ('acf_lags', '<u4'),
58 58 ('sea_to_atmos', '<f4'),
59 59 ('sea_notch', '<u4'),
60 60 ('lh_sea', '<u4'),
61 61 ('hh_sea', '<u4'),
62 62 ('nbins_sea', '<u4'),
63 63 ('min_snr', '<f4'),
64 64 ('min_cc', '<f4'),
65 65 ('max_time_diff', '<f4')
66 66 ])
67 67
68 68 DATA_STRUCTURE = numpy.dtype([
69 69 ('range', '<u4'),
70 70 ('status', '<u4'),
71 71 ('zonal', '<f4'),
72 72 ('meridional', '<f4'),
73 73 ('vertical', '<f4'),
74 74 ('zonal_a', '<f4'),
75 75 ('meridional_a', '<f4'),
76 76 ('corrected_fading', '<f4'), # seconds
77 77 ('uncorrected_fading', '<f4'), # seconds
78 78 ('time_diff', '<f4'),
79 79 ('major_axis', '<f4'),
80 80 ('axial_ratio', '<f4'),
81 81 ('orientation', '<f4'),
82 82 ('sea_power', '<u4'),
83 83 ('sea_algorithm', '<u4')
84 84 ])
85 85
86 86 @MPDecorator
87 87 class BLTRParamReader(JRODataReader, ProcessingUnit):
88 88 '''
89 89 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
90 90 '''
91 91
92 92 ext = '.sswma'
93 93
94 def __init__(self, **kwargs):
94 def __init__(self):
95 95
96 ProcessingUnit.__init__(self, **kwargs)
96 ProcessingUnit.__init__(self)
97 97
98 98 self.dataOut = Parameters()
99 99 self.counter_records = 0
100 100 self.flagNoMoreFiles = 0
101 101 self.isConfig = False
102 102 self.filename = None
103 103
104 104 def setup(self,
105 105 path=None,
106 106 startDate=None,
107 107 endDate=None,
108 108 ext=None,
109 109 startTime=datetime.time(0, 0, 0),
110 110 endTime=datetime.time(23, 59, 59),
111 111 timezone=0,
112 112 status_value=0,
113 113 **kwargs):
114 114
115 115 self.path = path
116 116 self.startDate = startDate
117 117 self.endDate = endDate
118 118 self.startTime = startTime
119 119 self.endTime = endTime
120 120 self.status_value = status_value
121 121 self.datatime = datetime.datetime(1900,1,1)
122 122
123 123 if self.path is None:
124 124 raise ValueError("The path is not valid")
125 125
126 126 if ext is None:
127 127 ext = self.ext
128 128
129 129 self.search_files(self.path, startDate, endDate, ext)
130 130 self.timezone = timezone
131 131 self.fileIndex = 0
132 132
133 133 if not self.fileList:
134 134 raise Warning("There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' " % (
135 135 path))
136 136
137 137 self.setNextFile()
138 138
139 139 def search_files(self, path, startDate, endDate, ext):
140 140 '''
141 141 Searching for BLTR rawdata file in path
142 142 Creating a list of file to proces included in [startDate,endDate]
143 143
144 144 Input:
145 145 path - Path to find BLTR rawdata files
146 146 startDate - Select file from this date
147 147 enDate - Select file until this date
148 148 ext - Extension of the file to read
149 149 '''
150 150
151 151 log.success('Searching files in {} '.format(path), 'BLTRParamReader')
152 152 foldercounter = 0
153 153 fileList0 = glob.glob1(path, "*%s" % ext)
154 154 fileList0.sort()
155 155
156 156 self.fileList = []
157 157 self.dateFileList = []
158 158
159 159 for thisFile in fileList0:
160 160 year = thisFile[-14:-10]
161 161 if not isNumber(year):
162 162 continue
163 163
164 164 month = thisFile[-10:-8]
165 165 if not isNumber(month):
166 166 continue
167 167
168 168 day = thisFile[-8:-6]
169 169 if not isNumber(day):
170 170 continue
171 171
172 172 year, month, day = int(year), int(month), int(day)
173 173 dateFile = datetime.date(year, month, day)
174 174
175 175 if (startDate > dateFile) or (endDate < dateFile):
176 176 continue
177 177
178 178 self.fileList.append(thisFile)
179 179 self.dateFileList.append(dateFile)
180 180
181 181 return
182 182
183 183 def setNextFile(self):
184 184
185 185 file_id = self.fileIndex
186 186
187 187 if file_id == len(self.fileList):
188 188 self.flagNoMoreFiles = 1
189 189 return 0
190 190
191 191 log.success('Opening {}'.format(self.fileList[file_id]), 'BLTRParamReader')
192 192 filename = os.path.join(self.path, self.fileList[file_id])
193 193
194 194 dirname, name = os.path.split(filename)
195 195 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
196 196 self.siteFile = name.split('.')[0]
197 197 if self.filename is not None:
198 198 self.fp.close()
199 199 self.filename = filename
200 200 self.fp = open(self.filename, 'rb')
201 201 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
202 202 self.nrecords = self.header_file['nrec'][0]
203 203 self.sizeOfFile = os.path.getsize(self.filename)
204 204 self.counter_records = 0
205 205 self.flagIsNewFile = 0
206 206 self.fileIndex += 1
207 207
208 208 return 1
209 209
210 210 def readNextBlock(self):
211 211
212 212 while True:
213 213 if self.counter_records == self.nrecords:
214 214 self.flagIsNewFile = 1
215 215 if not self.setNextFile():
216 216 return 0
217 217
218 218 self.readBlock()
219 219
220 220 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
221 221 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
222 222 log.warning(
223 223 'Reading Record No. {}/{} -> {} [Skipping]'.format(
224 224 self.counter_records,
225 225 self.nrecords,
226 226 self.datatime.ctime()),
227 227 'BLTRParamReader')
228 228 continue
229 229 break
230 230
231 231 log.log('Reading Record No. {}/{} -> {}'.format(
232 232 self.counter_records,
233 233 self.nrecords,
234 234 self.datatime.ctime()), 'BLTRParamReader')
235 235
236 236 return 1
237 237
238 238 def readBlock(self):
239 239
240 240 pointer = self.fp.tell()
241 241 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
242 242 self.nchannels = int(header_rec['nchan'][0] / 2)
243 243 self.kchan = header_rec['nrxs'][0]
244 244 self.nmodes = header_rec['nmodes'][0]
245 245 self.nranges = header_rec['nranges'][0]
246 246 self.fp.seek(pointer)
247 247 self.height = numpy.empty((self.nmodes, self.nranges))
248 self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
248 self.snr = numpy.empty((self.nmodes, int(self.nchannels), self.nranges))
249 249 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
250 250 self.flagDiscontinuousBlock = 0
251 251
252 252 for mode in range(self.nmodes):
253 253 self.readHeader()
254 254 data = self.readData()
255 255 self.height[mode] = (data[0] - self.correction) / 1000.
256 256 self.buffer[mode] = data[1]
257 257 self.snr[mode] = data[2]
258 258
259 259 self.counter_records = self.counter_records + self.nmodes
260 260
261 261 return
262 262
263 263 def readHeader(self):
264 264 '''
265 265 RecordHeader of BLTR rawdata file
266 266 '''
267 267
268 268 header_structure = numpy.dtype(
269 269 REC_HEADER_STRUCTURE.descr + [
270 ('antenna_coord', 'f4', (2, self.nchannels)),
271 ('rx_gains', 'u4', (self.nchannels,)),
272 ('rx_analysis', 'u4', (self.nchannels,))
270 ('antenna_coord', 'f4', (2, int(self.nchannels))),
271 ('rx_gains', 'u4', (int(self.nchannels),)),
272 ('rx_analysis', 'u4', (int(self.nchannels),))
273 273 ]
274 274 )
275 275
276 276 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
277 277 self.lat = self.header_rec['lat'][0]
278 278 self.lon = self.header_rec['lon'][0]
279 279 self.delta = self.header_rec['delta_r'][0]
280 280 self.correction = self.header_rec['dmode_rngcorr'][0]
281 281 self.imode = self.header_rec['dmode_index'][0]
282 282 self.antenna = self.header_rec['antenna_coord']
283 283 self.rx_gains = self.header_rec['rx_gains']
284 284 self.time = self.header_rec['time'][0]
285 285 dt = datetime.datetime.utcfromtimestamp(self.time)
286 286 if dt.date()>self.datatime.date():
287 287 self.flagDiscontinuousBlock = 1
288 288 self.datatime = dt
289 289
290 290 def readData(self):
291 291 '''
292 292 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
293 293
294 294 Input:
295 295 status_value - Array data is set to NAN for values that are not equal to status_value
296 296
297 297 '''
298 self.nchannels = int(self.nchannels)
298 299
299 300 data_structure = numpy.dtype(
300 301 DATA_STRUCTURE.descr + [
301 302 ('rx_saturation', 'u4', (self.nchannels,)),
302 303 ('chan_offset', 'u4', (2 * self.nchannels,)),
303 304 ('rx_amp', 'u4', (self.nchannels,)),
304 305 ('rx_snr', 'f4', (self.nchannels,)),
305 306 ('cross_snr', 'f4', (self.kchan,)),
306 307 ('sea_power_relative', 'f4', (self.kchan,))]
307 308 )
308 309
309 310 data = numpy.fromfile(self.fp, data_structure, self.nranges)
310 311
311 312 height = data['range']
312 313 winds = numpy.array(
313 314 (data['zonal'], data['meridional'], data['vertical']))
314 315 snr = data['rx_snr'].T
315 316
316 317 winds[numpy.where(winds == -9999.)] = numpy.nan
317 318 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
318 319 snr[numpy.where(snr == -9999.)] = numpy.nan
319 320 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
320 321 snr = numpy.power(10, snr / 10)
321 322
322 323 return height, winds, snr
323 324
324 325 def set_output(self):
325 326 '''
326 327 Storing data from databuffer to dataOut object
327 328 '''
328 329
329 330 self.dataOut.data_SNR = self.snr
330 331 self.dataOut.height = self.height
331 332 self.dataOut.data = self.buffer
332 333 self.dataOut.utctimeInit = self.time
333 334 self.dataOut.utctime = self.dataOut.utctimeInit
334 335 self.dataOut.useLocalTime = False
335 336 self.dataOut.paramInterval = 157
336 337 self.dataOut.timezone = self.timezone
337 338 self.dataOut.site = self.siteFile
338 339 self.dataOut.nrecords = self.nrecords / self.nmodes
339 340 self.dataOut.sizeOfFile = self.sizeOfFile
340 341 self.dataOut.lat = self.lat
341 342 self.dataOut.lon = self.lon
342 343 self.dataOut.channelList = list(range(self.nchannels))
343 344 self.dataOut.kchan = self.kchan
344 345 self.dataOut.delta = self.delta
345 346 self.dataOut.correction = self.correction
346 347 self.dataOut.nmodes = self.nmodes
347 348 self.dataOut.imode = self.imode
348 349 self.dataOut.antenna = self.antenna
349 350 self.dataOut.rx_gains = self.rx_gains
350 351 self.dataOut.flagNoData = False
351 352 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
352 353
353 354 def getData(self):
354 355 '''
355 356 Storing data from databuffer to dataOut object
356 357 '''
357 358 if self.flagNoMoreFiles:
358 359 self.dataOut.flagNoData = True
359 360 self.dataOut.error = (1, 'No More files to read')
360 361
361 362 if not self.readNextBlock():
362 363 self.dataOut.flagNoData = True
363 364 return 0
364 365
365 366 self.set_output()
366 367
367 return 1 No newline at end of file
368 return 1
369 No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,401 +1,402
1 1 '''
2 2 Created on Oct 24, 2016
3 3
4 4 @author: roj- LouVD
5 5 '''
6 6
7 7 import numpy
8 8 import copy
9 9 import datetime
10 10 import time
11 11 from time import gmtime
12 12
13 13 from numpy import transpose
14 14
15 15 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
16 16 from schainpy.model.data.jrodata import Parameters
17 17
18 18 @MPDecorator
19 19 class BLTRParametersProc(ProcessingUnit):
20 20 '''
21 21 Processing unit for BLTR parameters data (winds)
22 22
23 23 Inputs:
24 24 self.dataOut.nmodes - Number of operation modes
25 25 self.dataOut.nchannels - Number of channels
26 26 self.dataOut.nranges - Number of ranges
27 27
28 28 self.dataOut.data_SNR - SNR array
29 29 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
30 30 self.dataOut.height - Height array (km)
31 31 self.dataOut.time - Time array (seconds)
32 32
33 33 self.dataOut.fileIndex -Index of the file currently read
34 34 self.dataOut.lat - Latitude coordinate of BLTR location
35 35
36 36 self.dataOut.doy - Experiment doy (number of the day in the current year)
37 37 self.dataOut.month - Experiment month
38 38 self.dataOut.day - Experiment day
39 39 self.dataOut.year - Experiment year
40 40 '''
41 41
42 42 def __init__(self):
43 43 '''
44 44 Inputs: None
45 45 '''
46 46 ProcessingUnit.__init__(self)
47 47 self.dataOut = Parameters()
48 48
49 49 def setup(self, mode):
50 50 '''
51 51 '''
52 52 self.dataOut.mode = mode
53 53
54 54 def run(self, mode, snr_threshold=None):
55 55 '''
56 56 Inputs:
57 57 mode = High resolution (0) or Low resolution (1) data
58 58 snr_threshold = snr filter value
59 59 '''
60 60
61 61 if not self.isConfig:
62 62 self.setup(mode)
63 63 self.isConfig = True
64 64
65 65 if self.dataIn.type == 'Parameters':
66 66 self.dataOut.copy(self.dataIn)
67 67
68 68 self.dataOut.data_param = self.dataOut.data[mode]
69 69 self.dataOut.heightList = self.dataOut.height[0]
70 70 self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
71 71
72 72 if snr_threshold is not None:
73 73 SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
74 74 SNRavgdB = 10*numpy.log10(SNRavg)
75 75 for i in range(3):
76 76 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
77 77
78 78 # TODO
79 @MPDecorator
79 80 class OutliersFilter(Operation):
80 81
81 def __init__(self, **kwargs):
82 def __init__(self):
82 83 '''
83 84 '''
84 Operation.__init__(self, **kwargs)
85 Operation.__init__(self)
85 86
86 87 def run(self, svalue2, method, factor, filter, npoints=9):
87 88 '''
88 89 Inputs:
89 90 svalue - string to select array velocity
90 91 svalue2 - string to choose axis filtering
91 92 method - 0 for SMOOTH or 1 for MEDIAN
92 93 factor - number used to set threshold
93 94 filter - 1 for data filtering using the standard deviation criteria else 0
94 95 npoints - number of points for mask filter
95 96 '''
96 97
97 98 print(' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor))
98 99
99 100
100 101 yaxis = self.dataOut.heightList
101 102 xaxis = numpy.array([[self.dataOut.utctime]])
102 103
103 104 # Zonal
104 105 value_temp = self.dataOut.data_output[0]
105 106
106 107 # Zonal
107 108 value_temp = self.dataOut.data_output[1]
108 109
109 110 # Vertical
110 111 value_temp = numpy.transpose(self.dataOut.data_output[2])
111 112
112 113 htemp = yaxis
113 114 std = value_temp
114 115 for h in range(len(htemp)):
115 116 nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
116 117 minvalid = npoints
117 118
118 119 #only if valid values greater than the minimum required (10%)
119 120 if nvalues_valid > minvalid:
120 121
121 122 if method == 0:
122 123 #SMOOTH
123 124 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
124 125
125 126
126 127 if method == 1:
127 128 #MEDIAN
128 129 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
129 130
130 131 dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
131 132
132 133 threshold = dw*factor
133 134 value_temp[numpy.where(w > threshold),h] = numpy.nan
134 135 value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
135 136
136 137
137 138 #At the end
138 139 if svalue2 == 'inHeight':
139 140 value_temp = numpy.transpose(value_temp)
140 141 output_array[:,m] = value_temp
141 142
142 143 if svalue == 'zonal':
143 144 self.dataOut.data_output[0] = output_array
144 145
145 146 elif svalue == 'meridional':
146 147 self.dataOut.data_output[1] = output_array
147 148
148 149 elif svalue == 'vertical':
149 150 self.dataOut.data_output[2] = output_array
150 151
151 152 return self.dataOut.data_output
152 153
153 154
154 155 def Median(self,input,width):
155 156 '''
156 157 Inputs:
157 158 input - Velocity array
158 159 width - Number of points for mask filter
159 160
160 161 '''
161 162
162 163 if numpy.mod(width,2) == 1:
163 164 pc = int((width - 1) / 2)
164 165 cont = 0
165 166 output = []
166 167
167 168 for i in range(len(input)):
168 169 if i >= pc and i < len(input) - pc:
169 170 new2 = input[i-pc:i+pc+1]
170 171 temp = numpy.where(numpy.isfinite(new2))
171 172 new = new2[temp]
172 173 value = numpy.median(new)
173 174 output.append(value)
174 175
175 176 output = numpy.array(output)
176 177 output = numpy.hstack((input[0:pc],output))
177 178 output = numpy.hstack((output,input[-pc:len(input)]))
178 179
179 180 return output
180 181
181 182 def Smooth(self,input,width,edge_truncate = None):
182 183 '''
183 184 Inputs:
184 185 input - Velocity array
185 186 width - Number of points for mask filter
186 187 edge_truncate - 1 for truncate the convolution product else
187 188
188 189 '''
189 190
190 191 if numpy.mod(width,2) == 0:
191 192 real_width = width + 1
192 193 nzeros = width / 2
193 194 else:
194 195 real_width = width
195 196 nzeros = (width - 1) / 2
196 197
197 198 half_width = int(real_width)/2
198 199 length = len(input)
199 200
200 201 gate = numpy.ones(real_width,dtype='float')
201 202 norm_of_gate = numpy.sum(gate)
202 203
203 204 nan_process = 0
204 205 nan_id = numpy.where(numpy.isnan(input))
205 206 if len(nan_id[0]) > 0:
206 207 nan_process = 1
207 208 pb = numpy.zeros(len(input))
208 209 pb[nan_id] = 1.
209 210 input[nan_id] = 0.
210 211
211 212 if edge_truncate == True:
212 213 output = numpy.convolve(input/norm_of_gate,gate,mode='same')
213 214 elif edge_truncate == False or edge_truncate == None:
214 215 output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
215 216 output = numpy.hstack((input[0:half_width],output))
216 217 output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
217 218
218 219 if nan_process:
219 220 pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
220 221 pb = numpy.hstack((numpy.zeros(half_width),pb))
221 222 pb = numpy.hstack((pb,numpy.zeros(half_width)))
222 223 output[numpy.where(pb > 0.9999)] = numpy.nan
223 224 input[nan_id] = numpy.nan
224 225 return output
225 226
226 227 def Average(self,aver=0,nhaver=1):
227 228 '''
228 229 Inputs:
229 230 aver - Indicates the time period over which is averaged or consensus data
230 231 nhaver - Indicates the decimation factor in heights
231 232
232 233 '''
233 234 nhpoints = 48
234 235
235 236 lat_piura = -5.17
236 237 lat_huancayo = -12.04
237 238 lat_porcuya = -5.8
238 239
239 240 if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
240 241 hcm = 3.
241 242 if self.dataOut.year == 2003 :
242 243 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
243 244 nhpoints = 12
244 245
245 246 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
246 247 hcm = 3.
247 248 if self.dataOut.year == 2003 :
248 249 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
249 250 nhpoints = 12
250 251
251 252
252 253 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
253 254 hcm = 5.#2
254 255
255 256 pdata = 0.2
256 257 taver = [1,2,3,4,6,8,12,24]
257 258 t0 = 0
258 259 tf = 24
259 260 ntime =(tf-t0)/taver[aver]
260 261 ti = numpy.arange(ntime)
261 262 tf = numpy.arange(ntime) + taver[aver]
262 263
263 264
264 265 old_height = self.dataOut.heightList
265 266
266 267 if nhaver > 1:
267 268 num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
268 269 deltha = 0.05*nhaver
269 270 minhvalid = pdata*nhaver
270 271 for im in range(self.dataOut.nmodes):
271 272 new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
272 273
273 274
274 275 data_fHeigths_List = []
275 276 data_fZonal_List = []
276 277 data_fMeridional_List = []
277 278 data_fVertical_List = []
278 279 startDTList = []
279 280
280 281
281 282 for i in range(ntime):
282 283 height = old_height
283 284
284 285 start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
285 286 stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
286 287
287 288
288 289 limit_sec1 = time.mktime(start.timetuple())
289 290 limit_sec2 = time.mktime(stop.timetuple())
290 291
291 292 t1 = numpy.where(self.f_timesec >= limit_sec1)
292 293 t2 = numpy.where(self.f_timesec < limit_sec2)
293 294 time_select = []
294 295 for val_sec in t1[0]:
295 296 if val_sec in t2[0]:
296 297 time_select.append(val_sec)
297 298
298 299
299 300 time_select = numpy.array(time_select,dtype = 'int')
300 301 minvalid = numpy.ceil(pdata*nhpoints)
301 302
302 303 zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
303 304 mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
304 305 ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
305 306
306 307 if nhaver > 1:
307 308 new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
308 309 new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
309 310 new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
310 311
311 312 if len(time_select) > minvalid:
312 313 time_average = self.f_timesec[time_select]
313 314
314 315 for im in range(self.dataOut.nmodes):
315 316
316 317 for ih in range(self.dataOut.nranges):
317 318 if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
318 319 zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
319 320
320 321 if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
321 322 mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
322 323
323 324 if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
324 325 ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
325 326
326 327 if nhaver > 1:
327 328 for ih in range(num_hei):
328 329 hvalid = numpy.arange(nhaver) + nhaver*ih
329 330
330 331 if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
331 332 new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
332 333
333 334 if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
334 335 new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
335 336
336 337 if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
337 338 new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
338 339 if nhaver > 1:
339 340 zon_aver = new_zon_aver
340 341 mer_aver = new_mer_aver
341 342 ver_aver = new_ver_aver
342 343 height = new_height
343 344
344 345
345 346 tstart = time_average[0]
346 347 tend = time_average[-1]
347 348 startTime = time.gmtime(tstart)
348 349
349 350 year = startTime.tm_year
350 351 month = startTime.tm_mon
351 352 day = startTime.tm_mday
352 353 hour = startTime.tm_hour
353 354 minute = startTime.tm_min
354 355 second = startTime.tm_sec
355 356
356 357 startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
357 358
358 359
359 360 o_height = numpy.array([])
360 361 o_zon_aver = numpy.array([])
361 362 o_mer_aver = numpy.array([])
362 363 o_ver_aver = numpy.array([])
363 364 if self.dataOut.nmodes > 1:
364 365 for im in range(self.dataOut.nmodes):
365 366
366 367 if im == 0:
367 368 h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
368 369 else:
369 370 h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
370 371
371 372
372 373 ht = h_select[0]
373 374
374 375 o_height = numpy.hstack((o_height,height[im,ht]))
375 376 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
376 377 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
377 378 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
378 379
379 380 data_fHeigths_List.append(o_height)
380 381 data_fZonal_List.append(o_zon_aver)
381 382 data_fMeridional_List.append(o_mer_aver)
382 383 data_fVertical_List.append(o_ver_aver)
383 384
384 385
385 386 else:
386 387 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
387 388 ht = h_select[0]
388 389 o_height = numpy.hstack((o_height,height[im,ht]))
389 390 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
390 391 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
391 392 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
392 393
393 394 data_fHeigths_List.append(o_height)
394 395 data_fZonal_List.append(o_zon_aver)
395 396 data_fMeridional_List.append(o_mer_aver)
396 397 data_fVertical_List.append(o_ver_aver)
397 398
398 399
399 400 return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
400 401
401 402
Show file before | Show file after | Hide comments
@@ -1,951 +1,952
1 1 import itertools
2 2
3 3 import numpy
4 4
5 5 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
6 6 from schainpy.model.data.jrodata import Spectra
7 7 from schainpy.model.data.jrodata import hildebrand_sekhon
8 8 from schainpy.utils import log
9 9
10 10 @MPDecorator
11 11 class SpectraProc(ProcessingUnit):
12 12
13 13
14 14 def __init__(self):
15 15
16 16 ProcessingUnit.__init__(self)
17 17
18 18 self.buffer = None
19 19 self.firstdatatime = None
20 20 self.profIndex = 0
21 21 self.dataOut = Spectra()
22 22 self.id_min = None
23 23 self.id_max = None
24 24 self.setupReq = False #Agregar a todas las unidades de proc
25 25
26 26 def __updateSpecFromVoltage(self):
27 27
28 28 self.dataOut.timeZone = self.dataIn.timeZone
29 29 self.dataOut.dstFlag = self.dataIn.dstFlag
30 30 self.dataOut.errorCount = self.dataIn.errorCount
31 31 self.dataOut.useLocalTime = self.dataIn.useLocalTime
32 32 try:
33 33 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
34 34 except:
35 35 pass
36 36 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
37 37 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
38 38 self.dataOut.channelList = self.dataIn.channelList
39 39 self.dataOut.heightList = self.dataIn.heightList
40 40 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
41 41
42 42 self.dataOut.nBaud = self.dataIn.nBaud
43 43 self.dataOut.nCode = self.dataIn.nCode
44 44 self.dataOut.code = self.dataIn.code
45 45 self.dataOut.nProfiles = self.dataOut.nFFTPoints
46 46
47 47 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
48 48 self.dataOut.utctime = self.firstdatatime
49 49 # asumo q la data esta decodificada
50 50 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
51 51 # asumo q la data esta sin flip
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
53 53 self.dataOut.flagShiftFFT = False
54 54
55 55 self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 self.dataOut.nIncohInt = 1
57 57
58 58 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59
60 60 self.dataOut.frequency = self.dataIn.frequency
61 61 self.dataOut.realtime = self.dataIn.realtime
62 62
63 63 self.dataOut.azimuth = self.dataIn.azimuth
64 64 self.dataOut.zenith = self.dataIn.zenith
65 65
66 66 self.dataOut.beam.codeList = self.dataIn.beam.codeList
67 67 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
68 68 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
69 69
70 70 def __getFft(self):
71 71 """
72 72 Convierte valores de Voltaje a Spectra
73 73
74 74 Affected:
75 75 self.dataOut.data_spc
76 76 self.dataOut.data_cspc
77 77 self.dataOut.data_dc
78 78 self.dataOut.heightList
79 79 self.profIndex
80 80 self.buffer
81 81 self.dataOut.flagNoData
82 82 """
83 83 fft_volt = numpy.fft.fft(
84 84 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
85 85 fft_volt = fft_volt.astype(numpy.dtype('complex'))
86 86 dc = fft_volt[:, 0, :]
87 87
88 88 # calculo de self-spectra
89 89 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
90 90 spc = fft_volt * numpy.conjugate(fft_volt)
91 91 spc = spc.real
92 92
93 93 blocksize = 0
94 94 blocksize += dc.size
95 95 blocksize += spc.size
96 96
97 97 cspc = None
98 98 pairIndex = 0
99 99 if self.dataOut.pairsList != None:
100 100 # calculo de cross-spectra
101 101 cspc = numpy.zeros(
102 102 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
103 103 for pair in self.dataOut.pairsList:
104 104 if pair[0] not in self.dataOut.channelList:
105 105 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
106 106 str(pair), str(self.dataOut.channelList)))
107 107 if pair[1] not in self.dataOut.channelList:
108 108 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
109 109 str(pair), str(self.dataOut.channelList)))
110 110
111 111 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
112 112 numpy.conjugate(fft_volt[pair[1], :, :])
113 113 pairIndex += 1
114 114 blocksize += cspc.size
115 115
116 116 self.dataOut.data_spc = spc
117 117 self.dataOut.data_cspc = cspc
118 118 self.dataOut.data_dc = dc
119 119 self.dataOut.blockSize = blocksize
120 120 self.dataOut.flagShiftFFT = True
121 121
122 122 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None, shift_fft=False):
123 123
124 124 if self.dataIn.type == "Spectra":
125 125 self.dataOut.copy(self.dataIn)
126 126 if shift_fft:
127 127 #desplaza a la derecha en el eje 2 determinadas posiciones
128 128 shift = int(self.dataOut.nFFTPoints/2)
129 129 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
130 130
131 131 if self.dataOut.data_cspc is not None:
132 132 #desplaza a la derecha en el eje 2 determinadas posiciones
133 133 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
134 134
135 135 return True
136 136
137 137 if self.dataIn.type == "Voltage":
138 138
139 self.dataOut.flagNoData = True
140
139 141 if nFFTPoints == None:
140 142 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
141 143
142 144 if nProfiles == None:
143 145 nProfiles = nFFTPoints
144 146
145 147 if ippFactor == None:
146 148 ippFactor = 1
147 149
148 150 self.dataOut.ippFactor = ippFactor
149 151
150 152 self.dataOut.nFFTPoints = nFFTPoints
151 153 self.dataOut.pairsList = pairsList
152 154
153 155 if self.buffer is None:
154 156 self.buffer = numpy.zeros((self.dataIn.nChannels,
155 157 nProfiles,
156 158 self.dataIn.nHeights),
157 159 dtype='complex')
158 160
159 161 if self.dataIn.flagDataAsBlock:
160 162 # data dimension: [nChannels, nProfiles, nSamples]
161 163 nVoltProfiles = self.dataIn.data.shape[1]
162 164 # nVoltProfiles = self.dataIn.nProfiles
163 165
164 166 if nVoltProfiles == nProfiles:
165 167 self.buffer = self.dataIn.data.copy()
166 168 self.profIndex = nVoltProfiles
167 169
168 170 elif nVoltProfiles < nProfiles:
169 171
170 172 if self.profIndex == 0:
171 173 self.id_min = 0
172 174 self.id_max = nVoltProfiles
173 175
174 176 self.buffer[:, self.id_min:self.id_max,
175 177 :] = self.dataIn.data
176 178 self.profIndex += nVoltProfiles
177 179 self.id_min += nVoltProfiles
178 180 self.id_max += nVoltProfiles
179 181 else:
180 182 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
181 183 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
182 184 self.dataOut.flagNoData = True
183 185 return 0
184 186 else:
185 187 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
186 188 self.profIndex += 1
187 189
188 190 if self.firstdatatime == None:
189 191 self.firstdatatime = self.dataIn.utctime
190 192
191 193 if self.profIndex == nProfiles:
192 194 self.__updateSpecFromVoltage()
193 195 self.__getFft()
194 196
195 197 self.dataOut.flagNoData = False
196 198 self.firstdatatime = None
197 199 self.profIndex = 0
198 200
199 201 return True
200 202
201 203 raise ValueError("The type of input object '%s' is not valid" % (
202 204 self.dataIn.type))
203 205
204 206 def __selectPairs(self, pairsList):
205 207
206 208 if not pairsList:
207 209 return
208 210
209 211 pairs = []
210 212 pairsIndex = []
211 213
212 214 for pair in pairsList:
213 215 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
214 216 continue
215 217 pairs.append(pair)
216 218 pairsIndex.append(pairs.index(pair))
217 219
218 220 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
219 221 self.dataOut.pairsList = pairs
220 222
221 223 return
222 224
223 225 def __selectPairsByChannel(self, channelList=None):
224 226
225 227 if channelList == None:
226 228 return
227 229
228 230 pairsIndexListSelected = []
229 231 for pairIndex in self.dataOut.pairsIndexList:
230 232 # First pair
231 233 if self.dataOut.pairsList[pairIndex][0] not in channelList:
232 234 continue
233 235 # Second pair
234 236 if self.dataOut.pairsList[pairIndex][1] not in channelList:
235 237 continue
236 238
237 239 pairsIndexListSelected.append(pairIndex)
238 240
239 241 if not pairsIndexListSelected:
240 242 self.dataOut.data_cspc = None
241 243 self.dataOut.pairsList = []
242 244 return
243 245
244 246 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
245 247 self.dataOut.pairsList = [self.dataOut.pairsList[i]
246 248 for i in pairsIndexListSelected]
247 249
248 250 return
249 251
250 252 def selectChannels(self, channelList):
251 253
252 254 channelIndexList = []
253 255
254 256 for channel in channelList:
255 257 if channel not in self.dataOut.channelList:
256 258 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" % (
257 259 channel, str(self.dataOut.channelList)))
258 260
259 261 index = self.dataOut.channelList.index(channel)
260 262 channelIndexList.append(index)
261 263
262 264 self.selectChannelsByIndex(channelIndexList)
263 265
264 266 def selectChannelsByIndex(self, channelIndexList):
265 267 """
266 268 Selecciona un bloque de datos en base a canales segun el channelIndexList
267 269
268 270 Input:
269 271 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
270 272
271 273 Affected:
272 274 self.dataOut.data_spc
273 275 self.dataOut.channelIndexList
274 276 self.dataOut.nChannels
275 277
276 278 Return:
277 279 None
278 280 """
279 281
280 282 for channelIndex in channelIndexList:
281 283 if channelIndex not in self.dataOut.channelIndexList:
282 284 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " % (
283 285 channelIndex, self.dataOut.channelIndexList))
284 286
285 287 # nChannels = len(channelIndexList)
286 288
287 289 data_spc = self.dataOut.data_spc[channelIndexList, :]
288 290 data_dc = self.dataOut.data_dc[channelIndexList, :]
289 291
290 292 self.dataOut.data_spc = data_spc
291 293 self.dataOut.data_dc = data_dc
292 294
293 295 self.dataOut.channelList = [
294 296 self.dataOut.channelList[i] for i in channelIndexList]
295 297 # self.dataOut.nChannels = nChannels
296 298
297 299 self.__selectPairsByChannel(self.dataOut.channelList)
298 300
299 301 return 1
300 302
301 303 def selectHeights(self, minHei, maxHei):
302 304 """
303 305 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
304 306 minHei <= height <= maxHei
305 307
306 308 Input:
307 309 minHei : valor minimo de altura a considerar
308 310 maxHei : valor maximo de altura a considerar
309 311
310 312 Affected:
311 313 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
312 314
313 315 Return:
314 316 1 si el metodo se ejecuto con exito caso contrario devuelve 0
315 317 """
316 318
317 319 if (minHei > maxHei):
318 320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (
319 321 minHei, maxHei))
320 322
321 323 if (minHei < self.dataOut.heightList[0]):
322 324 minHei = self.dataOut.heightList[0]
323 325
324 326 if (maxHei > self.dataOut.heightList[-1]):
325 327 maxHei = self.dataOut.heightList[-1]
326 328
327 329 minIndex = 0
328 330 maxIndex = 0
329 331 heights = self.dataOut.heightList
330 332
331 333 inda = numpy.where(heights >= minHei)
332 334 indb = numpy.where(heights <= maxHei)
333 335
334 336 try:
335 337 minIndex = inda[0][0]
336 338 except:
337 339 minIndex = 0
338 340
339 341 try:
340 342 maxIndex = indb[0][-1]
341 343 except:
342 344 maxIndex = len(heights)
343 345
344 346 self.selectHeightsByIndex(minIndex, maxIndex)
345 347
346 348 return 1
347 349
348 350 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
349 351 newheis = numpy.where(
350 352 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
351 353
352 354 if hei_ref != None:
353 355 newheis = numpy.where(self.dataOut.heightList > hei_ref)
354 356
355 357 minIndex = min(newheis[0])
356 358 maxIndex = max(newheis[0])
357 359 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
358 360 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
359 361
360 362 # determina indices
361 363 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
362 364 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
363 365 avg_dB = 10 * \
364 366 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
365 367 beacon_dB = numpy.sort(avg_dB)[-nheis:]
366 368 beacon_heiIndexList = []
367 369 for val in avg_dB.tolist():
368 370 if val >= beacon_dB[0]:
369 371 beacon_heiIndexList.append(avg_dB.tolist().index(val))
370 372
371 373 #data_spc = data_spc[:,:,beacon_heiIndexList]
372 374 data_cspc = None
373 375 if self.dataOut.data_cspc is not None:
374 376 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
375 377 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
376 378
377 379 data_dc = None
378 380 if self.dataOut.data_dc is not None:
379 381 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
380 382 #data_dc = data_dc[:,beacon_heiIndexList]
381 383
382 384 self.dataOut.data_spc = data_spc
383 385 self.dataOut.data_cspc = data_cspc
384 386 self.dataOut.data_dc = data_dc
385 387 self.dataOut.heightList = heightList
386 388 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
387 389
388 390 return 1
389 391
390 392 def selectHeightsByIndex(self, minIndex, maxIndex):
391 393 """
392 394 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
393 395 minIndex <= index <= maxIndex
394 396
395 397 Input:
396 398 minIndex : valor de indice minimo de altura a considerar
397 399 maxIndex : valor de indice maximo de altura a considerar
398 400
399 401 Affected:
400 402 self.dataOut.data_spc
401 403 self.dataOut.data_cspc
402 404 self.dataOut.data_dc
403 405 self.dataOut.heightList
404 406
405 407 Return:
406 408 1 si el metodo se ejecuto con exito caso contrario devuelve 0
407 409 """
408 410
409 411 if (minIndex < 0) or (minIndex > maxIndex):
410 412 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
411 413 minIndex, maxIndex))
412 414
413 415 if (maxIndex >= self.dataOut.nHeights):
414 416 maxIndex = self.dataOut.nHeights - 1
415 417
416 418 # Spectra
417 419 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
418 420
419 421 data_cspc = None
420 422 if self.dataOut.data_cspc is not None:
421 423 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
422 424
423 425 data_dc = None
424 426 if self.dataOut.data_dc is not None:
425 427 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
426 428
427 429 self.dataOut.data_spc = data_spc
428 430 self.dataOut.data_cspc = data_cspc
429 431 self.dataOut.data_dc = data_dc
430 432
431 433 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
432 434
433 435 return 1
434 436
435 437 def removeDC(self, mode=2):
436 438 jspectra = self.dataOut.data_spc
437 439 jcspectra = self.dataOut.data_cspc
438 440
439 441 num_chan = jspectra.shape[0]
440 442 num_hei = jspectra.shape[2]
441 443
442 444 if jcspectra is not None:
443 445 jcspectraExist = True
444 446 num_pairs = jcspectra.shape[0]
445 447 else:
446 448 jcspectraExist = False
447 449
448 450 freq_dc = int(jspectra.shape[1] / 2)
449 451 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
450 452 ind_vel = ind_vel.astype(int)
451 453
452 454 if ind_vel[0] < 0:
453 455 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
454 456
455 457 if mode == 1:
456 458 jspectra[:, freq_dc, :] = (
457 459 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
458 460
459 461 if jcspectraExist:
460 462 jcspectra[:, freq_dc, :] = (
461 463 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
462 464
463 465 if mode == 2:
464 466
465 467 vel = numpy.array([-2, -1, 1, 2])
466 468 xx = numpy.zeros([4, 4])
467 469
468 470 for fil in range(4):
469 471 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
470 472
471 473 xx_inv = numpy.linalg.inv(xx)
472 474 xx_aux = xx_inv[0, :]
473 475
474 476 for ich in range(num_chan):
475 477 yy = jspectra[ich, ind_vel, :]
476 478 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
477 479
478 480 junkid = jspectra[ich, freq_dc, :] <= 0
479 481 cjunkid = sum(junkid)
480 482
481 483 if cjunkid.any():
482 484 jspectra[ich, freq_dc, junkid.nonzero()] = (
483 485 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
484 486
485 487 if jcspectraExist:
486 488 for ip in range(num_pairs):
487 489 yy = jcspectra[ip, ind_vel, :]
488 490 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
489 491
490 492 self.dataOut.data_spc = jspectra
491 493 self.dataOut.data_cspc = jcspectra
492 494
493 495 return 1
494 496
495 497 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
496 498
497 499 jspectra = self.dataOut.data_spc
498 500 jcspectra = self.dataOut.data_cspc
499 501 jnoise = self.dataOut.getNoise()
500 502 num_incoh = self.dataOut.nIncohInt
501 503
502 504 num_channel = jspectra.shape[0]
503 505 num_prof = jspectra.shape[1]
504 506 num_hei = jspectra.shape[2]
505 507
506 508 # hei_interf
507 509 if hei_interf is None:
508 510 count_hei = num_hei / 2 # Como es entero no importa
509 511 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
510 512 hei_interf = numpy.asarray(hei_interf)[0]
511 513 # nhei_interf
512 514 if (nhei_interf == None):
513 515 nhei_interf = 5
514 516 if (nhei_interf < 1):
515 517 nhei_interf = 1
516 518 if (nhei_interf > count_hei):
517 519 nhei_interf = count_hei
518 520 if (offhei_interf == None):
519 521 offhei_interf = 0
520 522
521 523 ind_hei = list(range(num_hei))
522 524 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
523 525 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
524 526 mask_prof = numpy.asarray(list(range(num_prof)))
525 527 num_mask_prof = mask_prof.size
526 528 comp_mask_prof = [0, num_prof / 2]
527 529
528 530 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
529 531 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
530 532 jnoise = numpy.nan
531 533 noise_exist = jnoise[0] < numpy.Inf
532 534
533 535 # Subrutina de Remocion de la Interferencia
534 536 for ich in range(num_channel):
535 537 # Se ordena los espectros segun su potencia (menor a mayor)
536 538 power = jspectra[ich, mask_prof, :]
537 539 power = power[:, hei_interf]
538 540 power = power.sum(axis=0)
539 541 psort = power.ravel().argsort()
540 542
541 543 # Se estima la interferencia promedio en los Espectros de Potencia empleando
542 544 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
543 545 offhei_interf, nhei_interf + offhei_interf))]]]
544 546
545 547 if noise_exist:
546 548 # tmp_noise = jnoise[ich] / num_prof
547 549 tmp_noise = jnoise[ich]
548 550 junkspc_interf = junkspc_interf - tmp_noise
549 551 #junkspc_interf[:,comp_mask_prof] = 0
550 552
551 553 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
552 554 jspc_interf = jspc_interf.transpose()
553 555 # Calculando el espectro de interferencia promedio
554 556 noiseid = numpy.where(
555 557 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
556 558 noiseid = noiseid[0]
557 559 cnoiseid = noiseid.size
558 560 interfid = numpy.where(
559 561 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
560 562 interfid = interfid[0]
561 563 cinterfid = interfid.size
562 564
563 565 if (cnoiseid > 0):
564 566 jspc_interf[noiseid] = 0
565 567
566 568 # Expandiendo los perfiles a limpiar
567 569 if (cinterfid > 0):
568 570 new_interfid = (
569 571 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
570 572 new_interfid = numpy.asarray(new_interfid)
571 573 new_interfid = {x for x in new_interfid}
572 574 new_interfid = numpy.array(list(new_interfid))
573 575 new_cinterfid = new_interfid.size
574 576 else:
575 577 new_cinterfid = 0
576 578
577 579 for ip in range(new_cinterfid):
578 580 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
579 581 jspc_interf[new_interfid[ip]
580 582 ] = junkspc_interf[ind[nhei_interf / 2], new_interfid[ip]]
581 583
582 584 jspectra[ich, :, ind_hei] = jspectra[ich, :,
583 585 ind_hei] - jspc_interf # Corregir indices
584 586
585 587 # Removiendo la interferencia del punto de mayor interferencia
586 588 ListAux = jspc_interf[mask_prof].tolist()
587 589 maxid = ListAux.index(max(ListAux))
588 590
589 591 if cinterfid > 0:
590 592 for ip in range(cinterfid * (interf == 2) - 1):
591 593 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
592 594 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
593 595 cind = len(ind)
594 596
595 597 if (cind > 0):
596 598 jspectra[ich, interfid[ip], ind] = tmp_noise * \
597 599 (1 + (numpy.random.uniform(cind) - 0.5) /
598 600 numpy.sqrt(num_incoh))
599 601
600 602 ind = numpy.array([-2, -1, 1, 2])
601 603 xx = numpy.zeros([4, 4])
602 604
603 605 for id1 in range(4):
604 606 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
605 607
606 608 xx_inv = numpy.linalg.inv(xx)
607 609 xx = xx_inv[:, 0]
608 610 ind = (ind + maxid + num_mask_prof) % num_mask_prof
609 611 yy = jspectra[ich, mask_prof[ind], :]
610 612 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
611 613 yy.transpose(), xx)
612 614
613 615 indAux = (jspectra[ich, :, :] < tmp_noise *
614 616 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
615 617 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
616 618 (1 - 1 / numpy.sqrt(num_incoh))
617 619
618 620 # Remocion de Interferencia en el Cross Spectra
619 621 if jcspectra is None:
620 622 return jspectra, jcspectra
621 623 num_pairs = jcspectra.size / (num_prof * num_hei)
622 624 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
623 625
624 626 for ip in range(num_pairs):
625 627
626 628 #-------------------------------------------
627 629
628 630 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
629 631 cspower = cspower[:, hei_interf]
630 632 cspower = cspower.sum(axis=0)
631 633
632 634 cspsort = cspower.ravel().argsort()
633 635 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
634 636 offhei_interf, nhei_interf + offhei_interf))]]]
635 637 junkcspc_interf = junkcspc_interf.transpose()
636 638 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
637 639
638 640 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
639 641
640 642 median_real = numpy.median(numpy.real(
641 643 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof / 4))]], :]))
642 644 median_imag = numpy.median(numpy.imag(
643 645 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof / 4))]], :]))
644 646 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
645 647 median_real, median_imag)
646 648
647 649 for iprof in range(num_prof):
648 650 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
649 651 jcspc_interf[iprof] = junkcspc_interf[iprof,
650 652 ind[nhei_interf / 2]]
651 653
652 654 # Removiendo la Interferencia
653 655 jcspectra[ip, :, ind_hei] = jcspectra[ip,
654 656 :, ind_hei] - jcspc_interf
655 657
656 658 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
657 659 maxid = ListAux.index(max(ListAux))
658 660
659 661 ind = numpy.array([-2, -1, 1, 2])
660 662 xx = numpy.zeros([4, 4])
661 663
662 664 for id1 in range(4):
663 665 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
664 666
665 667 xx_inv = numpy.linalg.inv(xx)
666 668 xx = xx_inv[:, 0]
667 669
668 670 ind = (ind + maxid + num_mask_prof) % num_mask_prof
669 671 yy = jcspectra[ip, mask_prof[ind], :]
670 672 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
671 673
672 674 # Guardar Resultados
673 675 self.dataOut.data_spc = jspectra
674 676 self.dataOut.data_cspc = jcspectra
675 677
676 678 return 1
677 679
678 680 def setRadarFrequency(self, frequency=None):
679 681
680 682 if frequency != None:
681 683 self.dataOut.frequency = frequency
682 684
683 685 return 1
684 686
685 687 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
686 688 # validacion de rango
687 689 if minHei == None:
688 690 minHei = self.dataOut.heightList[0]
689 691
690 692 if maxHei == None:
691 693 maxHei = self.dataOut.heightList[-1]
692 694
693 695 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
694 696 print('minHei: %.2f is out of the heights range' % (minHei))
695 697 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
696 698 minHei = self.dataOut.heightList[0]
697 699
698 700 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
699 701 print('maxHei: %.2f is out of the heights range' % (maxHei))
700 702 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
701 703 maxHei = self.dataOut.heightList[-1]
702 704
703 705 # validacion de velocidades
704 706 velrange = self.dataOut.getVelRange(1)
705 707
706 708 if minVel == None:
707 709 minVel = velrange[0]
708 710
709 711 if maxVel == None:
710 712 maxVel = velrange[-1]
711 713
712 714 if (minVel < velrange[0]) or (minVel > maxVel):
713 715 print('minVel: %.2f is out of the velocity range' % (minVel))
714 716 print('minVel is setting to %.2f' % (velrange[0]))
715 717 minVel = velrange[0]
716 718
717 719 if (maxVel > velrange[-1]) or (maxVel < minVel):
718 720 print('maxVel: %.2f is out of the velocity range' % (maxVel))
719 721 print('maxVel is setting to %.2f' % (velrange[-1]))
720 722 maxVel = velrange[-1]
721 723
722 724 # seleccion de indices para rango
723 725 minIndex = 0
724 726 maxIndex = 0
725 727 heights = self.dataOut.heightList
726 728
727 729 inda = numpy.where(heights >= minHei)
728 730 indb = numpy.where(heights <= maxHei)
729 731
730 732 try:
731 733 minIndex = inda[0][0]
732 734 except:
733 735 minIndex = 0
734 736
735 737 try:
736 738 maxIndex = indb[0][-1]
737 739 except:
738 740 maxIndex = len(heights)
739 741
740 742 if (minIndex < 0) or (minIndex > maxIndex):
741 743 raise ValueError("some value in (%d,%d) is not valid" % (
742 744 minIndex, maxIndex))
743 745
744 746 if (maxIndex >= self.dataOut.nHeights):
745 747 maxIndex = self.dataOut.nHeights - 1
746 748
747 749 # seleccion de indices para velocidades
748 750 indminvel = numpy.where(velrange >= minVel)
749 751 indmaxvel = numpy.where(velrange <= maxVel)
750 752 try:
751 753 minIndexVel = indminvel[0][0]
752 754 except:
753 755 minIndexVel = 0
754 756
755 757 try:
756 758 maxIndexVel = indmaxvel[0][-1]
757 759 except:
758 760 maxIndexVel = len(velrange)
759 761
760 762 # seleccion del espectro
761 763 data_spc = self.dataOut.data_spc[:,
762 764 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
763 765 # estimacion de ruido
764 766 noise = numpy.zeros(self.dataOut.nChannels)
765 767
766 768 for channel in range(self.dataOut.nChannels):
767 769 daux = data_spc[channel, :, :]
768 770 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
769 771
770 772 self.dataOut.noise_estimation = noise.copy()
771 773
772 774 return 1
773 775
774 776
775 777 class IncohInt(Operation):
776 778
777 779 __profIndex = 0
778 780 __withOverapping = False
779 781
780 782 __byTime = False
781 783 __initime = None
782 784 __lastdatatime = None
783 785 __integrationtime = None
784 786
785 787 __buffer_spc = None
786 788 __buffer_cspc = None
787 789 __buffer_dc = None
788 790
789 791 __dataReady = False
790 792
791 793 __timeInterval = None
792 794
793 795 n = None
794 796
795 797 def __init__(self):
796 798
797 799 Operation.__init__(self)
798 800
799 801 def setup(self, n=None, timeInterval=None, overlapping=False):
800 802 """
801 803 Set the parameters of the integration class.
802 804
803 805 Inputs:
804 806
805 807 n : Number of coherent integrations
806 808 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
807 809 overlapping :
808 810
809 811 """
810 812
811 813 self.__initime = None
812 814 self.__lastdatatime = 0
813 815
814 816 self.__buffer_spc = 0
815 817 self.__buffer_cspc = 0
816 818 self.__buffer_dc = 0
817 819
818 820 self.__profIndex = 0
819 821 self.__dataReady = False
820 822 self.__byTime = False
821 823
822 824 if n is None and timeInterval is None:
823 825 raise ValueError("n or timeInterval should be specified ...")
824 826
825 827 if n is not None:
826 828 self.n = int(n)
827 829 else:
828 830
829 831 self.__integrationtime = int(timeInterval)
830 832 self.n = None
831 833 self.__byTime = True
832 834
833 835 def putData(self, data_spc, data_cspc, data_dc):
834 836 """
835 837 Add a profile to the __buffer_spc and increase in one the __profileIndex
836 838
837 839 """
838 840
839 841 self.__buffer_spc += data_spc
840 842
841 843 if data_cspc is None:
842 844 self.__buffer_cspc = None
843 845 else:
844 846 self.__buffer_cspc += data_cspc
845 847
846 848 if data_dc is None:
847 849 self.__buffer_dc = None
848 850 else:
849 851 self.__buffer_dc += data_dc
850 852
851 853 self.__profIndex += 1
852 854
853 855 return
854 856
855 857 def pushData(self):
856 858 """
857 859 Return the sum of the last profiles and the profiles used in the sum.
858 860
859 861 Affected:
860 862
861 863 self.__profileIndex
862 864
863 865 """
864 866
865 867 data_spc = self.__buffer_spc
866 868 data_cspc = self.__buffer_cspc
867 869 data_dc = self.__buffer_dc
868 870 n = self.__profIndex
869 871
870 872 self.__buffer_spc = 0
871 873 self.__buffer_cspc = 0
872 874 self.__buffer_dc = 0
873 875 self.__profIndex = 0
874 876
875 877 return data_spc, data_cspc, data_dc, n
876 878
877 879 def byProfiles(self, *args):
878 880
879 881 self.__dataReady = False
880 882 avgdata_spc = None
881 883 avgdata_cspc = None
882 884 avgdata_dc = None
883 885
884 886 self.putData(*args)
885 887
886 888 if self.__profIndex == self.n:
887 889
888 890 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
889 891 self.n = n
890 892 self.__dataReady = True
891 893
892 894 return avgdata_spc, avgdata_cspc, avgdata_dc
893 895
894 896 def byTime(self, datatime, *args):
895 897
896 898 self.__dataReady = False
897 899 avgdata_spc = None
898 900 avgdata_cspc = None
899 901 avgdata_dc = None
900 902
901 903 self.putData(*args)
902 904
903 905 if (datatime - self.__initime) >= self.__integrationtime:
904 906 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
905 907 self.n = n
906 908 self.__dataReady = True
907 909
908 910 return avgdata_spc, avgdata_cspc, avgdata_dc
909 911
910 912 def integrate(self, datatime, *args):
911 913
912 914 if self.__profIndex == 0:
913 915 self.__initime = datatime
914 916
915 917 if self.__byTime:
916 918 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
917 919 datatime, *args)
918 920 else:
919 921 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
920 922
921 923 if not self.__dataReady:
922 924 return None, None, None, None
923 925
924 926 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
925 927
926 928 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
927 929 if n == 1:
928 930 return
929 931
930 932 dataOut.flagNoData = True
931 933
932 934 if not self.isConfig:
933 935 self.setup(n, timeInterval, overlapping)
934 936 self.isConfig = True
935 937
936 938 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
937 939 dataOut.data_spc,
938 940 dataOut.data_cspc,
939 941 dataOut.data_dc)
940 942
941 943 if self.__dataReady:
942 944
943 945 dataOut.data_spc = avgdata_spc
944 946 dataOut.data_cspc = avgdata_cspc
945 dataOut.data_dc = avgdata_dc
946
947 dataOut.data_dc = avgdata_dc
947 948 dataOut.nIncohInt *= self.n
948 949 dataOut.utctime = avgdatatime
949 950 dataOut.flagNoData = False
950 951
951 952 return dataOut No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1330 +1,1329
1 1 import sys
2 2 import numpy
3 3 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
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10 @MPDecorator
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 # self.dataOut.copy(self.dataIn)
30 30
31 31 def __updateObjFromAmisrInput(self):
32 32
33 33 self.dataOut.timeZone = self.dataIn.timeZone
34 34 self.dataOut.dstFlag = self.dataIn.dstFlag
35 35 self.dataOut.errorCount = self.dataIn.errorCount
36 36 self.dataOut.useLocalTime = self.dataIn.useLocalTime
37 37
38 38 self.dataOut.flagNoData = self.dataIn.flagNoData
39 39 self.dataOut.data = self.dataIn.data
40 40 self.dataOut.utctime = self.dataIn.utctime
41 41 self.dataOut.channelList = self.dataIn.channelList
42 42 #self.dataOut.timeInterval = self.dataIn.timeInterval
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.nProfiles = self.dataIn.nProfiles
45 45
46 46 self.dataOut.nCohInt = self.dataIn.nCohInt
47 47 self.dataOut.ippSeconds = self.dataIn.ippSeconds
48 48 self.dataOut.frequency = self.dataIn.frequency
49 49
50 50 self.dataOut.azimuth = self.dataIn.azimuth
51 51 self.dataOut.zenith = self.dataIn.zenith
52 52
53 53 self.dataOut.beam.codeList = self.dataIn.beam.codeList
54 54 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
55 55 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
56 56 #
57 57 # pass#
58 58 #
59 59 # def init(self):
60 60 #
61 61 #
62 62 # if self.dataIn.type == 'AMISR':
63 63 # self.__updateObjFromAmisrInput()
64 64 #
65 65 # if self.dataIn.type == 'Voltage':
66 66 # self.dataOut.copy(self.dataIn)
67 67 # # No necesita copiar en cada init() los atributos de dataIn
68 68 # # la copia deberia hacerse por cada nuevo bloque de datos
69 69
70 70 def selectChannels(self, channelList):
71 71
72 72 channelIndexList = []
73 73
74 74 for channel in channelList:
75 75 if channel not in self.dataOut.channelList:
76 76 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
77 77
78 78 index = self.dataOut.channelList.index(channel)
79 79 channelIndexList.append(index)
80 80
81 81 self.selectChannelsByIndex(channelIndexList)
82 82
83 83 def selectChannelsByIndex(self, channelIndexList):
84 84 """
85 85 Selecciona un bloque de datos en base a canales segun el channelIndexList
86 86
87 87 Input:
88 88 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
89 89
90 90 Affected:
91 91 self.dataOut.data
92 92 self.dataOut.channelIndexList
93 93 self.dataOut.nChannels
94 94 self.dataOut.m_ProcessingHeader.totalSpectra
95 95 self.dataOut.systemHeaderObj.numChannels
96 96 self.dataOut.m_ProcessingHeader.blockSize
97 97
98 98 Return:
99 99 None
100 100 """
101 101
102 102 for channelIndex in channelIndexList:
103 103 if channelIndex not in self.dataOut.channelIndexList:
104 104 print(channelIndexList)
105 105 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
106 106
107 107 if self.dataOut.flagDataAsBlock:
108 108 """
109 109 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
110 110 """
111 111 data = self.dataOut.data[channelIndexList,:,:]
112 112 else:
113 113 data = self.dataOut.data[channelIndexList,:]
114 114
115 115 self.dataOut.data = data
116 116 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
117 117 # self.dataOut.nChannels = nChannels
118 118
119 119 return 1
120 120
121 121 def selectHeights(self, minHei=None, maxHei=None):
122 122 """
123 123 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
124 124 minHei <= height <= maxHei
125 125
126 126 Input:
127 127 minHei : valor minimo de altura a considerar
128 128 maxHei : valor maximo de altura a considerar
129 129
130 130 Affected:
131 131 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
132 132
133 133 Return:
134 134 1 si el metodo se ejecuto con exito caso contrario devuelve 0
135 135 """
136 136
137 137 if minHei == None:
138 138 minHei = self.dataOut.heightList[0]
139 139
140 140 if maxHei == None:
141 141 maxHei = self.dataOut.heightList[-1]
142 142
143 143 if (minHei < self.dataOut.heightList[0]):
144 144 minHei = self.dataOut.heightList[0]
145 145
146 146 if (maxHei > self.dataOut.heightList[-1]):
147 147 maxHei = self.dataOut.heightList[-1]
148 148
149 149 minIndex = 0
150 150 maxIndex = 0
151 151 heights = self.dataOut.heightList
152 152
153 153 inda = numpy.where(heights >= minHei)
154 154 indb = numpy.where(heights <= maxHei)
155 155
156 156 try:
157 157 minIndex = inda[0][0]
158 158 except:
159 159 minIndex = 0
160 160
161 161 try:
162 162 maxIndex = indb[0][-1]
163 163 except:
164 164 maxIndex = len(heights)
165 165
166 166 self.selectHeightsByIndex(minIndex, maxIndex)
167 167
168 168 return 1
169 169
170 170
171 171 def selectHeightsByIndex(self, minIndex, maxIndex):
172 172 """
173 173 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
174 174 minIndex <= index <= maxIndex
175 175
176 176 Input:
177 177 minIndex : valor de indice minimo de altura a considerar
178 178 maxIndex : valor de indice maximo de altura a considerar
179 179
180 180 Affected:
181 181 self.dataOut.data
182 182 self.dataOut.heightList
183 183
184 184 Return:
185 185 1 si el metodo se ejecuto con exito caso contrario devuelve 0
186 186 """
187 187
188 188 if (minIndex < 0) or (minIndex > maxIndex):
189 189 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
190 190
191 191 if (maxIndex >= self.dataOut.nHeights):
192 192 maxIndex = self.dataOut.nHeights
193 193
194 194 #voltage
195 195 if self.dataOut.flagDataAsBlock:
196 196 """
197 197 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
198 198 """
199 199 data = self.dataOut.data[:,:, minIndex:maxIndex]
200 200 else:
201 201 data = self.dataOut.data[:, minIndex:maxIndex]
202 202
203 203 # firstHeight = self.dataOut.heightList[minIndex]
204 204
205 205 self.dataOut.data = data
206 206 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
207 207
208 208 if self.dataOut.nHeights <= 1:
209 209 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
210 210
211 211 return 1
212 212
213 213
214 214 def filterByHeights(self, window):
215 215
216 216 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
217 217
218 218 if window == None:
219 219 window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
220 220
221 221 newdelta = deltaHeight * window
222 222 r = self.dataOut.nHeights % window
223 223 newheights = (self.dataOut.nHeights-r)/window
224 224
225 225 if newheights <= 1:
226 226 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window))
227 227
228 228 if self.dataOut.flagDataAsBlock:
229 229 """
230 230 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
231 231 """
232 buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
232 buffer = self.dataOut.data[:, :, 0:int(self.dataOut.nHeights-r)]
233 233 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
234 234 buffer = numpy.sum(buffer,3)
235 235
236 236 else:
237 buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
238 buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
237 buffer = self.dataOut.data[:,0:int(self.dataOut.nHeights-r)]
238 buffer = buffer.reshape(self.dataOut.nChannels,int(self.dataOut.nHeights/window),int(window))
239 239 buffer = numpy.sum(buffer,2)
240 240
241 241 self.dataOut.data = buffer
242 242 self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
243 243 self.dataOut.windowOfFilter = window
244 244
245 245 def setH0(self, h0, deltaHeight = None):
246 246
247 247 if not deltaHeight:
248 248 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
249 249
250 250 nHeights = self.dataOut.nHeights
251 251
252 252 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
253 253
254 254 self.dataOut.heightList = newHeiRange
255 255
256 256 def deFlip(self, channelList = []):
257 257
258 258 data = self.dataOut.data.copy()
259 259
260 260 if self.dataOut.flagDataAsBlock:
261 261 flip = self.flip
262 262 profileList = list(range(self.dataOut.nProfiles))
263 263
264 264 if not channelList:
265 265 for thisProfile in profileList:
266 266 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
267 267 flip *= -1.0
268 268 else:
269 269 for thisChannel in channelList:
270 270 if thisChannel not in self.dataOut.channelList:
271 271 continue
272 272
273 273 for thisProfile in profileList:
274 274 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
275 275 flip *= -1.0
276 276
277 277 self.flip = flip
278 278
279 279 else:
280 280 if not channelList:
281 281 data[:,:] = data[:,:]*self.flip
282 282 else:
283 283 for thisChannel in channelList:
284 284 if thisChannel not in self.dataOut.channelList:
285 285 continue
286 286
287 287 data[thisChannel,:] = data[thisChannel,:]*self.flip
288 288
289 289 self.flip *= -1.
290 290
291 291 self.dataOut.data = data
292 292
293 293 def setRadarFrequency(self, frequency=None):
294 294
295 295 if frequency != None:
296 296 self.dataOut.frequency = frequency
297 297
298 298 return 1
299 299
300 300 def interpolateHeights(self, topLim, botLim):
301 301 #69 al 72 para julia
302 302 #82-84 para meteoros
303 303 if len(numpy.shape(self.dataOut.data))==2:
304 304 sampInterp = (self.dataOut.data[:,botLim-1] + self.dataOut.data[:,topLim+1])/2
305 305 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
306 306 #self.dataOut.data[:,botLim:limSup+1] = sampInterp
307 307 self.dataOut.data[:,botLim:topLim+1] = sampInterp
308 308 else:
309 309 nHeights = self.dataOut.data.shape[2]
310 310 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
311 311 y = self.dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
312 312 f = interpolate.interp1d(x, y, axis = 2)
313 313 xnew = numpy.arange(botLim,topLim+1)
314 314 ynew = f(xnew)
315 315
316 316 self.dataOut.data[:,:,botLim:topLim+1] = ynew
317 317
318 318 # import collections
319 319
320 320 class CohInt(Operation):
321 321
322 322 isConfig = False
323 323 __profIndex = 0
324 324 __byTime = False
325 325 __initime = None
326 326 __lastdatatime = None
327 327 __integrationtime = None
328 328 __buffer = None
329 329 __bufferStride = []
330 330 __dataReady = False
331 331 __profIndexStride = 0
332 332 __dataToPutStride = False
333 333 n = None
334 334
335 335 def __init__(self, **kwargs):
336 336
337 337 Operation.__init__(self, **kwargs)
338 338
339 339 # self.isConfig = False
340 340
341 341 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
342 342 """
343 343 Set the parameters of the integration class.
344 344
345 345 Inputs:
346 346
347 347 n : Number of coherent integrations
348 348 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
349 349 overlapping :
350 350 """
351 351
352 352 self.__initime = None
353 353 self.__lastdatatime = 0
354 354 self.__buffer = None
355 355 self.__dataReady = False
356 356 self.byblock = byblock
357 357 self.stride = stride
358 358
359 359 if n == None and timeInterval == None:
360 360 raise ValueError("n or timeInterval should be specified ...")
361 361
362 362 if n != None:
363 363 self.n = n
364 364 self.__byTime = False
365 365 else:
366 366 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
367 367 self.n = 9999
368 368 self.__byTime = True
369 369
370 370 if overlapping:
371 371 self.__withOverlapping = True
372 372 self.__buffer = None
373 373 else:
374 374 self.__withOverlapping = False
375 375 self.__buffer = 0
376 376
377 377 self.__profIndex = 0
378 378
379 379 def putData(self, data):
380 380
381 381 """
382 382 Add a profile to the __buffer and increase in one the __profileIndex
383 383
384 384 """
385 385
386 386 if not self.__withOverlapping:
387 387 self.__buffer += data.copy()
388 388 self.__profIndex += 1
389 389 return
390 390
391 391 #Overlapping data
392 392 nChannels, nHeis = data.shape
393 393 data = numpy.reshape(data, (1, nChannels, nHeis))
394 394
395 395 #If the buffer is empty then it takes the data value
396 396 if self.__buffer is None:
397 397 self.__buffer = data
398 398 self.__profIndex += 1
399 399 return
400 400
401 401 #If the buffer length is lower than n then stakcing the data value
402 402 if self.__profIndex < self.n:
403 403 self.__buffer = numpy.vstack((self.__buffer, data))
404 404 self.__profIndex += 1
405 405 return
406 406
407 407 #If the buffer length is equal to n then replacing the last buffer value with the data value
408 408 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
409 409 self.__buffer[self.n-1] = data
410 410 self.__profIndex = self.n
411 411 return
412 412
413 413
414 414 def pushData(self):
415 415 """
416 416 Return the sum of the last profiles and the profiles used in the sum.
417 417
418 418 Affected:
419 419
420 420 self.__profileIndex
421 421
422 422 """
423 423
424 424 if not self.__withOverlapping:
425 425 data = self.__buffer
426 426 n = self.__profIndex
427 427
428 428 self.__buffer = 0
429 429 self.__profIndex = 0
430 430
431 431 return data, n
432 432
433 433 #Integration with Overlapping
434 434 data = numpy.sum(self.__buffer, axis=0)
435 435 # print data
436 436 # raise
437 437 n = self.__profIndex
438 438
439 439 return data, n
440 440
441 441 def byProfiles(self, data):
442 442
443 443 self.__dataReady = False
444 444 avgdata = None
445 445 # n = None
446 446 # print data
447 447 # raise
448 448 self.putData(data)
449 449
450 450 if self.__profIndex == self.n:
451 451 avgdata, n = self.pushData()
452 452 self.__dataReady = True
453 453
454 454 return avgdata
455 455
456 456 def byTime(self, data, datatime):
457 457
458 458 self.__dataReady = False
459 459 avgdata = None
460 460 n = None
461 461
462 462 self.putData(data)
463 463
464 464 if (datatime - self.__initime) >= self.__integrationtime:
465 465 avgdata, n = self.pushData()
466 466 self.n = n
467 467 self.__dataReady = True
468 468
469 469 return avgdata
470 470
471 471 def integrateByStride(self, data, datatime):
472 472 # print data
473 473 if self.__profIndex == 0:
474 474 self.__buffer = [[data.copy(), datatime]]
475 475 else:
476 476 self.__buffer.append([data.copy(),datatime])
477 477 self.__profIndex += 1
478 478 self.__dataReady = False
479 479
480 480 if self.__profIndex == self.n * self.stride :
481 481 self.__dataToPutStride = True
482 482 self.__profIndexStride = 0
483 483 self.__profIndex = 0
484 484 self.__bufferStride = []
485 485 for i in range(self.stride):
486 486 current = self.__buffer[i::self.stride]
487 487 data = numpy.sum([t[0] for t in current], axis=0)
488 488 avgdatatime = numpy.average([t[1] for t in current])
489 489 # print data
490 490 self.__bufferStride.append((data, avgdatatime))
491 491
492 492 if self.__dataToPutStride:
493 493 self.__dataReady = True
494 494 self.__profIndexStride += 1
495 495 if self.__profIndexStride == self.stride:
496 496 self.__dataToPutStride = False
497 497 # print self.__bufferStride[self.__profIndexStride - 1]
498 498 # raise
499 499 return self.__bufferStride[self.__profIndexStride - 1]
500 500
501 501
502 502 return None, None
503 503
504 504 def integrate(self, data, datatime=None):
505 505
506 506 if self.__initime == None:
507 507 self.__initime = datatime
508 508
509 509 if self.__byTime:
510 510 avgdata = self.byTime(data, datatime)
511 511 else:
512 512 avgdata = self.byProfiles(data)
513 513
514 514
515 515 self.__lastdatatime = datatime
516 516
517 517 if avgdata is None:
518 518 return None, None
519 519
520 520 avgdatatime = self.__initime
521 521
522 522 deltatime = datatime - self.__lastdatatime
523 523
524 524 if not self.__withOverlapping:
525 525 self.__initime = datatime
526 526 else:
527 527 self.__initime += deltatime
528 528
529 529 return avgdata, avgdatatime
530 530
531 531 def integrateByBlock(self, dataOut):
532 532
533 533 times = int(dataOut.data.shape[1]/self.n)
534 534 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
535 535
536 536 id_min = 0
537 537 id_max = self.n
538 538
539 539 for i in range(times):
540 540 junk = dataOut.data[:,id_min:id_max,:]
541 541 avgdata[:,i,:] = junk.sum(axis=1)
542 542 id_min += self.n
543 543 id_max += self.n
544 544
545 545 timeInterval = dataOut.ippSeconds*self.n
546 546 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
547 547 self.__dataReady = True
548 548 return avgdata, avgdatatime
549 549
550 550 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
551 551
552 552 if not self.isConfig:
553 553 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
554 554 self.isConfig = True
555 555
556 556 if dataOut.flagDataAsBlock:
557 557 """
558 558 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
559 559 """
560 560 avgdata, avgdatatime = self.integrateByBlock(dataOut)
561 561 dataOut.nProfiles /= self.n
562 562 else:
563 563 if stride is None:
564 564 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
565 565 else:
566 566 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
567 567
568 568
569 569 # dataOut.timeInterval *= n
570 570 dataOut.flagNoData = True
571 571
572 572 if self.__dataReady:
573 573 dataOut.data = avgdata
574 574 dataOut.nCohInt *= self.n
575 575 dataOut.utctime = avgdatatime
576 576 # print avgdata, avgdatatime
577 577 # raise
578 578 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
579 579 dataOut.flagNoData = False
580 580 return dataOut
581 581
582 582 class Decoder(Operation):
583 583
584 584 isConfig = False
585 585 __profIndex = 0
586 586
587 587 code = None
588 588
589 589 nCode = None
590 590 nBaud = None
591 591
592 592 def __init__(self, **kwargs):
593 593
594 594 Operation.__init__(self, **kwargs)
595 595
596 596 self.times = None
597 597 self.osamp = None
598 598 # self.__setValues = False
599 599 self.isConfig = False
600 600 self.setupReq = False
601 601 def setup(self, code, osamp, dataOut):
602 602
603 603 self.__profIndex = 0
604 604
605 605 self.code = code
606 606
607 607 self.nCode = len(code)
608 608 self.nBaud = len(code[0])
609 609
610 610 if (osamp != None) and (osamp >1):
611 611 self.osamp = osamp
612 612 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
613 613 self.nBaud = self.nBaud*self.osamp
614 614
615 615 self.__nChannels = dataOut.nChannels
616 616 self.__nProfiles = dataOut.nProfiles
617 617 self.__nHeis = dataOut.nHeights
618 618
619 619 if self.__nHeis < self.nBaud:
620 620 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
621 621
622 622 #Frequency
623 623 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
624 624
625 625 __codeBuffer[:,0:self.nBaud] = self.code
626 626
627 627 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
628 628
629 629 if dataOut.flagDataAsBlock:
630 630
631 631 self.ndatadec = self.__nHeis #- self.nBaud + 1
632 632
633 633 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
634 634
635 635 else:
636 636
637 637 #Time
638 638 self.ndatadec = self.__nHeis #- self.nBaud + 1
639 639
640 640 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
641 641
642 642 def __convolutionInFreq(self, data):
643 643
644 644 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
645 645
646 646 fft_data = numpy.fft.fft(data, axis=1)
647 647
648 648 conv = fft_data*fft_code
649 649
650 650 data = numpy.fft.ifft(conv,axis=1)
651 651
652 652 return data
653 653
654 654 def __convolutionInFreqOpt(self, data):
655 655
656 656 raise NotImplementedError
657 657
658 658 def __convolutionInTime(self, data):
659 659
660 660 code = self.code[self.__profIndex]
661 661 for i in range(self.__nChannels):
662 662 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
663 663
664 664 return self.datadecTime
665 665
666 666 def __convolutionByBlockInTime(self, data):
667 667
668 668 repetitions = self.__nProfiles / self.nCode
669 669
670 670 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
671 671 junk = junk.flatten()
672 672 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
673 673 profilesList = range(self.__nProfiles)
674 674
675 675 for i in range(self.__nChannels):
676 676 for j in profilesList:
677 677 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
678 678 return self.datadecTime
679 679
680 680 def __convolutionByBlockInFreq(self, data):
681 681
682 682 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
683 683
684 684
685 685 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
686 686
687 687 fft_data = numpy.fft.fft(data, axis=2)
688 688
689 689 conv = fft_data*fft_code
690 690
691 691 data = numpy.fft.ifft(conv,axis=2)
692 692
693 693 return data
694 694
695 695
696 696 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
697 697
698 698 if dataOut.flagDecodeData:
699 699 print("This data is already decoded, recoding again ...")
700 700
701 701 if not self.isConfig:
702 702
703 703 if code is None:
704 704 if dataOut.code is None:
705 705 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
706 706
707 707 code = dataOut.code
708 708 else:
709 709 code = numpy.array(code).reshape(nCode,nBaud)
710 710 self.setup(code, osamp, dataOut)
711 711
712 712 self.isConfig = True
713 713
714 714 if mode == 3:
715 715 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
716 716
717 717 if times != None:
718 718 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
719 719
720 720 if self.code is None:
721 721 print("Fail decoding: Code is not defined.")
722 722 return
723 723
724 724 self.__nProfiles = dataOut.nProfiles
725 725 datadec = None
726 726
727 727 if mode == 3:
728 728 mode = 0
729 729
730 730 if dataOut.flagDataAsBlock:
731 731 """
732 732 Decoding when data have been read as block,
733 733 """
734 734
735 735 if mode == 0:
736 736 datadec = self.__convolutionByBlockInTime(dataOut.data)
737 737 if mode == 1:
738 738 datadec = self.__convolutionByBlockInFreq(dataOut.data)
739 739 else:
740 740 """
741 741 Decoding when data have been read profile by profile
742 742 """
743 743 if mode == 0:
744 744 datadec = self.__convolutionInTime(dataOut.data)
745 745
746 746 if mode == 1:
747 747 datadec = self.__convolutionInFreq(dataOut.data)
748 748
749 749 if mode == 2:
750 750 datadec = self.__convolutionInFreqOpt(dataOut.data)
751 751
752 752 if datadec is None:
753 753 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
754 754
755 755 dataOut.code = self.code
756 756 dataOut.nCode = self.nCode
757 757 dataOut.nBaud = self.nBaud
758 758
759 759 dataOut.data = datadec
760 760
761 761 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
762 762
763 763 dataOut.flagDecodeData = True #asumo q la data esta decodificada
764 764
765 765 if self.__profIndex == self.nCode-1:
766 766 self.__profIndex = 0
767 767 return dataOut
768 768
769 769 self.__profIndex += 1
770 770
771 771 return dataOut
772 772 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
773 773
774 774
775 775 class ProfileConcat(Operation):
776 776
777 777 isConfig = False
778 778 buffer = None
779 779
780 780 def __init__(self, **kwargs):
781 781
782 782 Operation.__init__(self, **kwargs)
783 783 self.profileIndex = 0
784 784
785 785 def reset(self):
786 786 self.buffer = numpy.zeros_like(self.buffer)
787 787 self.start_index = 0
788 788 self.times = 1
789 789
790 790 def setup(self, data, m, n=1):
791 791 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
792 792 self.nHeights = data.shape[1]#.nHeights
793 793 self.start_index = 0
794 794 self.times = 1
795 795
796 796 def concat(self, data):
797 797
798 798 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
799 799 self.start_index = self.start_index + self.nHeights
800 800
801 801 def run(self, dataOut, m):
802
803 802 dataOut.flagNoData = True
804 803
805 804 if not self.isConfig:
806 805 self.setup(dataOut.data, m, 1)
807 806 self.isConfig = True
808 807
809 808 if dataOut.flagDataAsBlock:
810 809 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
811 810
812 811 else:
813 812 self.concat(dataOut.data)
814 813 self.times += 1
815 814 if self.times > m:
816 815 dataOut.data = self.buffer
817 816 self.reset()
818 817 dataOut.flagNoData = False
819 818 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
820 819 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
821 820 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
822 821 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
823 822 dataOut.ippSeconds *= m
824 823 return dataOut
825 824
826 825 class ProfileSelector(Operation):
827 826
828 827 profileIndex = None
829 828 # Tamanho total de los perfiles
830 829 nProfiles = None
831 830
832 831 def __init__(self, **kwargs):
833 832
834 833 Operation.__init__(self, **kwargs)
835 834 self.profileIndex = 0
836 835
837 836 def incProfileIndex(self):
838 837
839 838 self.profileIndex += 1
840 839
841 840 if self.profileIndex >= self.nProfiles:
842 841 self.profileIndex = 0
843 842
844 843 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
845 844
846 845 if profileIndex < minIndex:
847 846 return False
848 847
849 848 if profileIndex > maxIndex:
850 849 return False
851 850
852 851 return True
853 852
854 853 def isThisProfileInList(self, profileIndex, profileList):
855 854
856 855 if profileIndex not in profileList:
857 856 return False
858 857
859 858 return True
860 859
861 860 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
862 861
863 862 """
864 863 ProfileSelector:
865 864
866 865 Inputs:
867 866 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
868 867
869 868 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
870 869
871 870 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
872 871
873 872 """
874 873
875 874 if rangeList is not None:
876 875 if type(rangeList[0]) not in (tuple, list):
877 876 rangeList = [rangeList]
878 877
879 878 dataOut.flagNoData = True
880 879
881 880 if dataOut.flagDataAsBlock:
882 881 """
883 882 data dimension = [nChannels, nProfiles, nHeis]
884 883 """
885 884 if profileList != None:
886 885 dataOut.data = dataOut.data[:,profileList,:]
887 886
888 887 if profileRangeList != None:
889 888 minIndex = profileRangeList[0]
890 889 maxIndex = profileRangeList[1]
891 890 profileList = list(range(minIndex, maxIndex+1))
892 891
893 892 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
894 893
895 894 if rangeList != None:
896 895
897 896 profileList = []
898 897
899 898 for thisRange in rangeList:
900 899 minIndex = thisRange[0]
901 900 maxIndex = thisRange[1]
902 901
903 902 profileList.extend(list(range(minIndex, maxIndex+1)))
904 903
905 904 dataOut.data = dataOut.data[:,profileList,:]
906 905
907 906 dataOut.nProfiles = len(profileList)
908 907 dataOut.profileIndex = dataOut.nProfiles - 1
909 908 dataOut.flagNoData = False
910 909
911 return True
910 return dataOut
912 911
913 912 """
914 913 data dimension = [nChannels, nHeis]
915 914 """
916 915
917 916 if profileList != None:
918 917
919 918 if self.isThisProfileInList(dataOut.profileIndex, profileList):
920 919
921 920 self.nProfiles = len(profileList)
922 921 dataOut.nProfiles = self.nProfiles
923 922 dataOut.profileIndex = self.profileIndex
924 923 dataOut.flagNoData = False
925 924
926 925 self.incProfileIndex()
927 return True
926 return dataOut
928 927
929 928 if profileRangeList != None:
930 929
931 930 minIndex = profileRangeList[0]
932 931 maxIndex = profileRangeList[1]
933 932
934 933 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
935 934
936 935 self.nProfiles = maxIndex - minIndex + 1
937 936 dataOut.nProfiles = self.nProfiles
938 937 dataOut.profileIndex = self.profileIndex
939 938 dataOut.flagNoData = False
940 939
941 940 self.incProfileIndex()
942 return True
941 return dataOut
943 942
944 943 if rangeList != None:
945 944
946 945 nProfiles = 0
947 946
948 947 for thisRange in rangeList:
949 948 minIndex = thisRange[0]
950 949 maxIndex = thisRange[1]
951 950
952 951 nProfiles += maxIndex - minIndex + 1
953 952
954 953 for thisRange in rangeList:
955 954
956 955 minIndex = thisRange[0]
957 956 maxIndex = thisRange[1]
958 957
959 958 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
960 959
961 960 self.nProfiles = nProfiles
962 961 dataOut.nProfiles = self.nProfiles
963 962 dataOut.profileIndex = self.profileIndex
964 963 dataOut.flagNoData = False
965 964
966 965 self.incProfileIndex()
967 966
968 967 break
969 968
970 return True
969 return dataOut
971 970
972 971
973 972 if beam != None: #beam is only for AMISR data
974 973 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
975 974 dataOut.flagNoData = False
976 975 dataOut.profileIndex = self.profileIndex
977 976
978 977 self.incProfileIndex()
979 978
980 return True
979 return dataOut
981 980
982 981 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
983 982
984 983 #return False
985 984 return dataOut
986 985
987 986 class Reshaper(Operation):
988 987
989 988 def __init__(self, **kwargs):
990 989
991 990 Operation.__init__(self, **kwargs)
992 991
993 992 self.__buffer = None
994 993 self.__nitems = 0
995 994
996 995 def __appendProfile(self, dataOut, nTxs):
997 996
998 997 if self.__buffer is None:
999 998 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1000 999 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1001 1000
1002 1001 ini = dataOut.nHeights * self.__nitems
1003 1002 end = ini + dataOut.nHeights
1004 1003
1005 1004 self.__buffer[:, ini:end] = dataOut.data
1006 1005
1007 1006 self.__nitems += 1
1008 1007
1009 1008 return int(self.__nitems*nTxs)
1010 1009
1011 1010 def __getBuffer(self):
1012 1011
1013 1012 if self.__nitems == int(1./self.__nTxs):
1014 1013
1015 1014 self.__nitems = 0
1016 1015
1017 1016 return self.__buffer.copy()
1018 1017
1019 1018 return None
1020 1019
1021 1020 def __checkInputs(self, dataOut, shape, nTxs):
1022 1021
1023 1022 if shape is None and nTxs is None:
1024 1023 raise ValueError("Reshaper: shape of factor should be defined")
1025 1024
1026 1025 if nTxs:
1027 1026 if nTxs < 0:
1028 1027 raise ValueError("nTxs should be greater than 0")
1029 1028
1030 1029 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1031 1030 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1032 1031
1033 1032 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1034 1033
1035 1034 return shape, nTxs
1036 1035
1037 1036 if len(shape) != 2 and len(shape) != 3:
1038 1037 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1039 1038
1040 1039 if len(shape) == 2:
1041 1040 shape_tuple = [dataOut.nChannels]
1042 1041 shape_tuple.extend(shape)
1043 1042 else:
1044 1043 shape_tuple = list(shape)
1045 1044
1046 1045 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1047 1046
1048 1047 return shape_tuple, nTxs
1049 1048
1050 1049 def run(self, dataOut, shape=None, nTxs=None):
1051 1050
1052 1051 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1053 1052
1054 1053 dataOut.flagNoData = True
1055 1054 profileIndex = None
1056 1055
1057 1056 if dataOut.flagDataAsBlock:
1058 1057
1059 1058 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1060 1059 dataOut.flagNoData = False
1061 1060
1062 1061 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1063 1062
1064 1063 else:
1065 1064
1066 1065 if self.__nTxs < 1:
1067 1066
1068 1067 self.__appendProfile(dataOut, self.__nTxs)
1069 1068 new_data = self.__getBuffer()
1070 1069
1071 1070 if new_data is not None:
1072 1071 dataOut.data = new_data
1073 1072 dataOut.flagNoData = False
1074 1073
1075 1074 profileIndex = dataOut.profileIndex*nTxs
1076 1075
1077 1076 else:
1078 1077 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1079 1078
1080 1079 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1081 1080
1082 1081 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1083 1082
1084 1083 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1085 1084
1086 1085 dataOut.profileIndex = profileIndex
1087 1086
1088 1087 dataOut.ippSeconds /= self.__nTxs
1089 1088
1090 1089 return dataOut
1091 1090
1092 1091 class SplitProfiles(Operation):
1093 1092
1094 1093 def __init__(self, **kwargs):
1095 1094
1096 1095 Operation.__init__(self, **kwargs)
1097 1096
1098 1097 def run(self, dataOut, n):
1099 1098
1100 1099 dataOut.flagNoData = True
1101 1100 profileIndex = None
1102 1101
1103 1102 if dataOut.flagDataAsBlock:
1104 1103
1105 1104 #nchannels, nprofiles, nsamples
1106 1105 shape = dataOut.data.shape
1107 1106
1108 1107 if shape[2] % n != 0:
1109 1108 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1110 1109
1111 1110 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1112 1111
1113 1112 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1114 1113 dataOut.flagNoData = False
1115 1114
1116 1115 profileIndex = int(dataOut.nProfiles/n) - 1
1117 1116
1118 1117 else:
1119 1118
1120 1119 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1121 1120
1122 1121 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1123 1122
1124 1123 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1125 1124
1126 1125 dataOut.nProfiles = int(dataOut.nProfiles*n)
1127 1126
1128 1127 dataOut.profileIndex = profileIndex
1129 1128
1130 1129 dataOut.ippSeconds /= n
1131 1130
1132 1131 return dataOut
1133 1132
1134 1133 class CombineProfiles(Operation):
1135 1134 def __init__(self, **kwargs):
1136 1135
1137 1136 Operation.__init__(self, **kwargs)
1138 1137
1139 1138 self.__remData = None
1140 1139 self.__profileIndex = 0
1141 1140
1142 1141 def run(self, dataOut, n):
1143 1142
1144 1143 dataOut.flagNoData = True
1145 1144 profileIndex = None
1146 1145
1147 1146 if dataOut.flagDataAsBlock:
1148 1147
1149 1148 #nchannels, nprofiles, nsamples
1150 1149 shape = dataOut.data.shape
1151 1150 new_shape = shape[0], shape[1]/n, shape[2]*n
1152 1151
1153 1152 if shape[1] % n != 0:
1154 1153 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1155 1154
1156 1155 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1157 1156 dataOut.flagNoData = False
1158 1157
1159 1158 profileIndex = int(dataOut.nProfiles*n) - 1
1160 1159
1161 1160 else:
1162 1161
1163 1162 #nchannels, nsamples
1164 1163 if self.__remData is None:
1165 1164 newData = dataOut.data
1166 1165 else:
1167 1166 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1168 1167
1169 1168 self.__profileIndex += 1
1170 1169
1171 1170 if self.__profileIndex < n:
1172 1171 self.__remData = newData
1173 1172 #continue
1174 1173 return
1175 1174
1176 1175 self.__profileIndex = 0
1177 1176 self.__remData = None
1178 1177
1179 1178 dataOut.data = newData
1180 1179 dataOut.flagNoData = False
1181 1180
1182 1181 profileIndex = dataOut.profileIndex/n
1183 1182
1184 1183
1185 1184 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1186 1185
1187 1186 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1188 1187
1189 1188 dataOut.nProfiles = int(dataOut.nProfiles/n)
1190 1189
1191 1190 dataOut.profileIndex = profileIndex
1192 1191
1193 1192 dataOut.ippSeconds *= n
1194 1193
1195 1194 return dataOut
1196 1195 # import collections
1197 1196 # from scipy.stats import mode
1198 1197 #
1199 1198 # class Synchronize(Operation):
1200 1199 #
1201 1200 # isConfig = False
1202 1201 # __profIndex = 0
1203 1202 #
1204 1203 # def __init__(self, **kwargs):
1205 1204 #
1206 1205 # Operation.__init__(self, **kwargs)
1207 1206 # # self.isConfig = False
1208 1207 # self.__powBuffer = None
1209 1208 # self.__startIndex = 0
1210 1209 # self.__pulseFound = False
1211 1210 #
1212 1211 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1213 1212 #
1214 1213 # #Read data
1215 1214 #
1216 1215 # powerdB = dataOut.getPower(channel = channel)
1217 1216 # noisedB = dataOut.getNoise(channel = channel)[0]
1218 1217 #
1219 1218 # self.__powBuffer.extend(powerdB.flatten())
1220 1219 #
1221 1220 # dataArray = numpy.array(self.__powBuffer)
1222 1221 #
1223 1222 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1224 1223 #
1225 1224 # maxValue = numpy.nanmax(filteredPower)
1226 1225 #
1227 1226 # if maxValue < noisedB + 10:
1228 1227 # #No se encuentra ningun pulso de transmision
1229 1228 # return None
1230 1229 #
1231 1230 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1232 1231 #
1233 1232 # if len(maxValuesIndex) < 2:
1234 1233 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1235 1234 # return None
1236 1235 #
1237 1236 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1238 1237 #
1239 1238 # #Seleccionar solo valores con un espaciamiento de nSamples
1240 1239 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1241 1240 #
1242 1241 # if len(pulseIndex) < 2:
1243 1242 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1244 1243 # return None
1245 1244 #
1246 1245 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1247 1246 #
1248 1247 # #remover senales que se distancien menos de 10 unidades o muestras
1249 1248 # #(No deberian existir IPP menor a 10 unidades)
1250 1249 #
1251 1250 # realIndex = numpy.where(spacing > 10 )[0]
1252 1251 #
1253 1252 # if len(realIndex) < 2:
1254 1253 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1255 1254 # return None
1256 1255 #
1257 1256 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1258 1257 # realPulseIndex = pulseIndex[realIndex]
1259 1258 #
1260 1259 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1261 1260 #
1262 1261 # print "IPP = %d samples" %period
1263 1262 #
1264 1263 # self.__newNSamples = dataOut.nHeights #int(period)
1265 1264 # self.__startIndex = int(realPulseIndex[0])
1266 1265 #
1267 1266 # return 1
1268 1267 #
1269 1268 #
1270 1269 # def setup(self, nSamples, nChannels, buffer_size = 4):
1271 1270 #
1272 1271 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1273 1272 # maxlen = buffer_size*nSamples)
1274 1273 #
1275 1274 # bufferList = []
1276 1275 #
1277 1276 # for i in range(nChannels):
1278 1277 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1279 1278 # maxlen = buffer_size*nSamples)
1280 1279 #
1281 1280 # bufferList.append(bufferByChannel)
1282 1281 #
1283 1282 # self.__nSamples = nSamples
1284 1283 # self.__nChannels = nChannels
1285 1284 # self.__bufferList = bufferList
1286 1285 #
1287 1286 # def run(self, dataOut, channel = 0):
1288 1287 #
1289 1288 # if not self.isConfig:
1290 1289 # nSamples = dataOut.nHeights
1291 1290 # nChannels = dataOut.nChannels
1292 1291 # self.setup(nSamples, nChannels)
1293 1292 # self.isConfig = True
1294 1293 #
1295 1294 # #Append new data to internal buffer
1296 1295 # for thisChannel in range(self.__nChannels):
1297 1296 # bufferByChannel = self.__bufferList[thisChannel]
1298 1297 # bufferByChannel.extend(dataOut.data[thisChannel])
1299 1298 #
1300 1299 # if self.__pulseFound:
1301 1300 # self.__startIndex -= self.__nSamples
1302 1301 #
1303 1302 # #Finding Tx Pulse
1304 1303 # if not self.__pulseFound:
1305 1304 # indexFound = self.__findTxPulse(dataOut, channel)
1306 1305 #
1307 1306 # if indexFound == None:
1308 1307 # dataOut.flagNoData = True
1309 1308 # return
1310 1309 #
1311 1310 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1312 1311 # self.__pulseFound = True
1313 1312 # self.__startIndex = indexFound
1314 1313 #
1315 1314 # #If pulse was found ...
1316 1315 # for thisChannel in range(self.__nChannels):
1317 1316 # bufferByChannel = self.__bufferList[thisChannel]
1318 1317 # #print self.__startIndex
1319 1318 # x = numpy.array(bufferByChannel)
1320 1319 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1321 1320 #
1322 1321 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1323 1322 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1324 1323 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1325 1324 #
1326 1325 # dataOut.data = self.__arrayBuffer
1327 1326 #
1328 1327 # self.__startIndex += self.__newNSamples
1329 1328 #
1330 1329 # return
General Comments 0
You need to be logged in to leave comments. Login now