##// END OF EJS Templates
schain
RSS
Adicionada operación Clean Rayleigh a Spectra Proc
joabAM -
r1385:83c0b1494611
parent child
Show More
Side by Side Unified
Context file:
r1385:83c0b1494611 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,101 +1,102
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plo Specra Heis data
6 6
7 7 """
8 8
9 9 import numpy
10 10
11 11 from schainpy.model.graphics.jroplot_base import Plot, plt
12
12 import matplotlib.pyplot as plt
13 13
14 14 class SpectraHeisPlot(Plot):
15 15
16 16 CODE = 'spc_heis'
17 17
18 18 def setup(self):
19 19
20 20 self.nplots = len(self.data.channels)
21 21 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
22 22 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
23 23 self.height = 2.6 * self.nrows
24 24 self.width = 3.5 * self.ncols
25 25 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.95, 'bottom': 0.08})
26 26 self.ylabel = 'Intensity [dB]'
27 27 self.xlabel = 'Frequency [KHz]'
28 28 self.colorbar = False
29 29
30 30 def update(self, dataOut):
31 31
32 32 data = {}
33 33 meta = {}
34 34 spc = 10*numpy.log10(dataOut.data_spc / dataOut.normFactor)
35 35 data['spc_heis'] = spc
36 36
37 37 return data, meta
38 38
39 39 def plot(self):
40 40
41 41 c = 3E8
42 deltaHeight = self.data.yrange[1] - self.data.yrange[0]
42 deltaHeight = self.data.yrange[1] - self.data.yrange[0] # yrange = heightList
43 43 x = numpy.arange(-1*len(self.data.yrange)/2., len(self.data.yrange)/2.)*(c/(2*deltaHeight*len(self.data.yrange)*1000))
44 44 self.y = self.data[-1]['spc_heis']
45 45 self.titles = []
46 46
47 Maintitle = "Range from %d km to %d km" %(int(self.data.yrange[0]),int(self.data.yrange[-1]))
47 48 for n, ax in enumerate(self.axes):
48 49 ychannel = self.y[n,:]
49 50 if ax.firsttime:
50 51 self.xmin = min(x) if self.xmin is None else self.xmin
51 52 self.xmax = max(x) if self.xmax is None else self.xmax
52 53 ax.plt = ax.plot(x, ychannel, lw=1, color='b')[0]
53 54 else:
54 55 ax.plt.set_data(x, ychannel)
55 56
56 57 self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
57
58 plt.suptitle(Maintitle)
58 59
59 60 class RTIHeisPlot(Plot):
60 61
61 62 CODE = 'rti_heis'
62 63
63 64 def setup(self):
64 65
65 66 self.xaxis = 'time'
66 67 self.ncols = 1
67 68 self.nrows = 1
68 69 self.nplots = 1
69 70 self.ylabel = 'Intensity [dB]'
70 71 self.xlabel = 'Time'
71 72 self.titles = ['RTI']
72 73 self.colorbar = False
73 74 self.height = 4
74 75 self.plots_adjust.update({'right': 0.85 })
75 76
76 77 def update(self, dataOut):
77 78
78 79 data = {}
79 80 meta = {}
80 81 spc = dataOut.data_spc / dataOut.normFactor
81 82 spc = 10*numpy.log10(numpy.average(spc, axis=1))
82 83 data['rti_heis'] = spc
83 84
84 85 return data, meta
85 86
86 87 def plot(self):
87 88
88 89 x = self.data.times
89 90 Y = self.data['rti_heis']
90 91
91 92 if self.axes[0].firsttime:
92 93 self.ymin = numpy.nanmin(Y) - 5 if self.ymin == None else self.ymin
93 94 self.ymax = numpy.nanmax(Y) + 5 if self.ymax == None else self.ymax
94 95 for ch in self.data.channels:
95 96 y = Y[ch]
96 97 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
97 98 plt.legend(bbox_to_anchor=(1.18, 1.0))
98 99 else:
99 100 for ch in self.data.channels:
100 101 y = Y[ch]
101 102 self.axes[0].lines[ch].set_data(x, y)
Show file before | Show file after | Hide comments
@@ -1,302 +1,302
1 1 '''
2 2 Created on Jul 9, 2014
3 3
4 4 @author: roj-idl71
5 5 '''
6 6 import os
7 7 import datetime
8 8 import numpy
9 9
10 10 from schainpy.model.graphics.jroplot_base import Plot, plt
11 11
12 12
13 13 class ScopePlot(Plot):
14 14
15 15 '''
16 16 Plot for Scope
17 17 '''
18 18
19 19 CODE = 'scope'
20 20 plot_type = 'scatter'
21 21
22 22 def setup(self):
23 23
24 24 self.xaxis = 'Range (Km)'
25 self.ncols = 1
26 self.nrows = 1
27 self.nplots = 1
25 self.nplots = len(self.data.channels)
26 self.nrows = int(numpy.ceil(self.nplots/2))
27 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
28 28 self.ylabel = 'Intensity [dB]'
29 self.titles = ['Scope']
29 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
30 30 self.colorbar = False
31 31 self.width = 6
32 32 self.height = 4
33 33
34 34 def update(self, dataOut):
35 35
36 36 data = {}
37 37 meta = {
38 38 'nProfiles': dataOut.nProfiles,
39 39 'flagDataAsBlock': dataOut.flagDataAsBlock,
40 40 'profileIndex': dataOut.profileIndex,
41 41 }
42 42 if self.CODE == 'scope':
43 43 data[self.CODE] = dataOut.data
44 44 elif self.CODE == 'pp_power':
45 45 data[self.CODE] = dataOut.dataPP_POWER
46 46 elif self.CODE == 'pp_signal':
47 47 data[self.CODE] = dataOut.dataPP_POW
48 48 elif self.CODE == 'pp_velocity':
49 49 data[self.CODE] = dataOut.dataPP_DOP
50 50 elif self.CODE == 'pp_specwidth':
51 51 data[self.CODE] = dataOut.dataPP_WIDTH
52 52
53 53 return data, meta
54 54
55 55 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
56 56
57 57 yreal = y[channelIndexList,:].real
58 58 yimag = y[channelIndexList,:].imag
59 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
59 Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
60 60 self.xlabel = "Range (Km)"
61 61 self.ylabel = "Intensity - IQ"
62 62
63 63 self.y = yreal
64 64 self.x = x
65 65
66 self.titles[0] = title
67
68 66 for i,ax in enumerate(self.axes):
69 67 title = "Channel %d" %(i)
70 68 if ax.firsttime:
71 69 self.xmin = min(x)
72 70 self.xmax = max(x)
73 71 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
74 72 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
75 73 else:
76 74 ax.plt_r.set_data(x, yreal[i,:])
77 75 ax.plt_i.set_data(x, yimag[i,:])
76 plt.suptitle(Maintitle)
78 77
79 78 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
80 79 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
81 80 yreal = y.real
82 81 yreal = 10*numpy.log10(yreal)
83 82 self.y = yreal
84 83 title = wintitle + " Power: %s" %(thisDatetime.strftime("%d-%b-%Y"))
85 84 self.xlabel = "Range (Km)"
86 85 self.ylabel = "Intensity [dB]"
87 86
88 87
89 88 self.titles[0] = title
90 89
91 90 for i,ax in enumerate(self.axes):
92 91 title = "Channel %d" %(i)
93 92 ychannel = yreal[i,:]
94 93
95 94 if ax.firsttime:
96 95 self.xmin = min(x)
97 96 self.xmax = max(x)
98 97 ax.plt_r = ax.plot(x, ychannel)[0]
99 98 else:
100 99 ax.plt_r.set_data(x, ychannel)
101 100
101
102 102 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
103 103
104 104
105 105 y = y[channelIndexList,:]
106 106 yreal = y.real
107 107 yreal = 10*numpy.log10(yreal)
108 108 self.y = yreal
109 109 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
110 110 self.xlabel = "Range (Km)"
111 111 self.ylabel = "Intensity"
112 112 self.xmin = min(x)
113 113 self.xmax = max(x)
114 114
115 115 self.titles[0] =title
116 116 for i,ax in enumerate(self.axes):
117 117 title = "Channel %d" %(i)
118 118
119 119 ychannel = yreal[i,:]
120 120
121 121 if ax.firsttime:
122 122 ax.plt_r = ax.plot(x, ychannel)[0]
123 123 else:
124 124 #pass
125 125 ax.plt_r.set_data(x, ychannel)
126 126
127 127 def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle):
128 128
129 129 x = x[channelIndexList,:]
130 130 yreal = y
131 131 self.y = yreal
132 132 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
133 133 self.xlabel = "Velocity (m/s)"
134 134 self.ylabel = "Range (Km)"
135 135 self.xmin = numpy.min(x)
136 136 self.xmax = numpy.max(x)
137 137 self.titles[0] =title
138 138 for i,ax in enumerate(self.axes):
139 139 title = "Channel %d" %(i)
140 140 xchannel = x[i,:]
141 141 if ax.firsttime:
142 142 ax.plt_r = ax.plot(xchannel, yreal)[0]
143 143 else:
144 144 #pass
145 145 ax.plt_r.set_data(xchannel, yreal)
146 146
147 147 def plot_weatherspecwidth(self, x, y, channelIndexList, thisDatetime, wintitle):
148 148
149 149 x = x[channelIndexList,:]
150 150 yreal = y
151 151 self.y = yreal
152 152 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
153 153 self.xlabel = "width "
154 154 self.ylabel = "Range (Km)"
155 155 self.xmin = numpy.min(x)
156 156 self.xmax = numpy.max(x)
157 157 self.titles[0] =title
158 158 for i,ax in enumerate(self.axes):
159 159 title = "Channel %d" %(i)
160 160 xchannel = x[i,:]
161 161 if ax.firsttime:
162 162 ax.plt_r = ax.plot(xchannel, yreal)[0]
163 163 else:
164 164 #pass
165 165 ax.plt_r.set_data(xchannel, yreal)
166 166
167 167 def plot(self):
168 168 if self.channels:
169 169 channels = self.channels
170 170 else:
171 171 channels = self.data.channels
172 172
173 173 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
174 174
175 175 scope = self.data[-1][self.CODE]
176 176
177 177 if self.data.flagDataAsBlock:
178 178
179 179 for i in range(self.data.nProfiles):
180 180
181 181 wintitle1 = " [Profile = %d] " %i
182 182 if self.CODE =="scope":
183 183 if self.type == "power":
184 184 self.plot_power(self.data.yrange,
185 185 scope[:,i,:],
186 186 channels,
187 187 thisDatetime,
188 188 wintitle1
189 189 )
190 190
191 191 if self.type == "iq":
192 192 self.plot_iq(self.data.yrange,
193 193 scope[:,i,:],
194 194 channels,
195 195 thisDatetime,
196 196 wintitle1
197 197 )
198 198 if self.CODE=="pp_power":
199 199 self.plot_weatherpower(self.data.yrange,
200 200 scope[:,i,:],
201 201 channels,
202 202 thisDatetime,
203 203 wintitle
204 204 )
205 205 if self.CODE=="pp_signal":
206 206 self.plot_weatherpower(self.data.yrange,
207 207 scope[:,i,:],
208 208 channels,
209 209 thisDatetime,
210 210 wintitle
211 211 )
212 212 if self.CODE=="pp_velocity":
213 213 self.plot_weathervelocity(scope[:,i,:],
214 214 self.data.yrange,
215 215 channels,
216 216 thisDatetime,
217 217 wintitle
218 218 )
219 219 if self.CODE=="pp_spcwidth":
220 220 self.plot_weatherspecwidth(scope[:,i,:],
221 221 self.data.yrange,
222 222 channels,
223 223 thisDatetime,
224 224 wintitle
225 225 )
226 226 else:
227 227 wintitle = " [Profile = %d] " %self.data.profileIndex
228 228 if self.CODE== "scope":
229 229 if self.type == "power":
230 230 self.plot_power(self.data.yrange,
231 231 scope,
232 232 channels,
233 233 thisDatetime,
234 234 wintitle
235 235 )
236 236
237 237 if self.type == "iq":
238 238 self.plot_iq(self.data.yrange,
239 239 scope,
240 240 channels,
241 241 thisDatetime,
242 242 wintitle
243 243 )
244 244 if self.CODE=="pp_power":
245 245 self.plot_weatherpower(self.data.yrange,
246 246 scope,
247 247 channels,
248 248 thisDatetime,
249 249 wintitle
250 250 )
251 251 if self.CODE=="pp_signal":
252 252 self.plot_weatherpower(self.data.yrange,
253 253 scope,
254 254 channels,
255 255 thisDatetime,
256 256 wintitle
257 257 )
258 258 if self.CODE=="pp_velocity":
259 259 self.plot_weathervelocity(scope,
260 260 self.data.yrange,
261 261 channels,
262 262 thisDatetime,
263 263 wintitle
264 264 )
265 265 if self.CODE=="pp_specwidth":
266 266 self.plot_weatherspecwidth(scope,
267 267 self.data.yrange,
268 268 channels,
269 269 thisDatetime,
270 270 wintitle
271 271 )
272 272
273 273
274 274 class PulsepairPowerPlot(ScopePlot):
275 275 '''
276 276 Plot for P= S+N
277 277 '''
278 278
279 279 CODE = 'pp_power'
280 280 plot_type = 'scatter'
281 281
282 282 class PulsepairVelocityPlot(ScopePlot):
283 283 '''
284 284 Plot for VELOCITY
285 285 '''
286 286 CODE = 'pp_velocity'
287 287 plot_type = 'scatter'
288 288
289 289 class PulsepairSpecwidthPlot(ScopePlot):
290 290 '''
291 291 Plot for WIDTH
292 292 '''
293 293 CODE = 'pp_specwidth'
294 294 plot_type = 'scatter'
295 295
296 296 class PulsepairSignalPlot(ScopePlot):
297 297 '''
298 298 Plot for S
299 299 '''
300 300
301 301 CODE = 'pp_signal'
302 302 plot_type = 'scatter'
Show file before | Show file after | Hide comments
@@ -1,1575 +1,1575
1 1 """
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 """
6 6 import os
7 7 import sys
8 8 import glob
9 9 import time
10 10 import numpy
11 11 import fnmatch
12 12 import inspect
13 13 import time
14 14 import datetime
15 15 import zmq
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, MPDecorator
18 18 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
19 19 from schainpy.model.data.jroheaderIO import get_dtype_index, get_numpy_dtype, get_procflag_dtype, get_dtype_width
20 20 from schainpy.utils import log
21 21 import schainpy.admin
22 22
23 23 LOCALTIME = True
24 24 DT_DIRECTIVES = {
25 25 '%Y': 4,
26 26 '%y': 2,
27 27 '%m': 2,
28 28 '%d': 2,
29 29 '%j': 3,
30 30 '%H': 2,
31 31 '%M': 2,
32 32 '%S': 2,
33 33 '%f': 6
34 34 }
35 35
36 36
37 37 def isNumber(cad):
38 38 """
39 39 Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero.
40 40
41 41 Excepciones:
42 42 Si un determinado string no puede ser convertido a numero
43 43 Input:
44 44 str, string al cual se le analiza para determinar si convertible a un numero o no
45 45
46 46 Return:
47 47 True : si el string es uno numerico
48 48 False : no es un string numerico
49 49 """
50 50 try:
51 51 float(cad)
52 52 return True
53 53 except:
54 54 return False
55 55
56 56
57 57 def isFileInEpoch(filename, startUTSeconds, endUTSeconds):
58 58 """
59 59 Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado.
60 60
61 61 Inputs:
62 62 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
63 63
64 64 startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en
65 65 segundos contados desde 01/01/1970.
66 66 endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en
67 67 segundos contados desde 01/01/1970.
68 68
69 69 Return:
70 70 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
71 71 fecha especificado, de lo contrario retorna False.
72 72
73 73 Excepciones:
74 74 Si el archivo no existe o no puede ser abierto
75 75 Si la cabecera no puede ser leida.
76 76
77 77 """
78 78 basicHeaderObj = BasicHeader(LOCALTIME)
79 79
80 80 try:
81 81 fp = open(filename, 'rb')
82 82 except IOError:
83 83 print("The file %s can't be opened" % (filename))
84 84 return 0
85 85
86 86 sts = basicHeaderObj.read(fp)
87 87 fp.close()
88 88
89 89 if not(sts):
90 90 print("Skipping the file %s because it has not a valid header" % (filename))
91 91 return 0
92 92
93 93 if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)):
94 94 return 0
95 95
96 96 return 1
97 97
98 98
99 99 def isTimeInRange(thisTime, startTime, endTime):
100 100 if endTime >= startTime:
101 101 if (thisTime < startTime) or (thisTime > endTime):
102 102 return 0
103 103 return 1
104 104 else:
105 105 if (thisTime < startTime) and (thisTime > endTime):
106 106 return 0
107 107 return 1
108 108
109 109
110 110 def isFileInTimeRange(filename, startDate, endDate, startTime, endTime):
111 111 """
112 112 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
113 113
114 114 Inputs:
115 115 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
116 116
117 117 startDate : fecha inicial del rango seleccionado en formato datetime.date
118 118
119 119 endDate : fecha final del rango seleccionado en formato datetime.date
120 120
121 121 startTime : tiempo inicial del rango seleccionado en formato datetime.time
122 122
123 123 endTime : tiempo final del rango seleccionado en formato datetime.time
124 124
125 125 Return:
126 126 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
127 127 fecha especificado, de lo contrario retorna False.
128 128
129 129 Excepciones:
130 130 Si el archivo no existe o no puede ser abierto
131 131 Si la cabecera no puede ser leida.
132 132
133 133 """
134 134
135 135 try:
136 136 fp = open(filename, 'rb')
137 137 except IOError:
138 138 print("The file %s can't be opened" % (filename))
139 139 return None
140 140
141 141 firstBasicHeaderObj = BasicHeader(LOCALTIME)
142 142 systemHeaderObj = SystemHeader()
143 143 radarControllerHeaderObj = RadarControllerHeader()
144 144 processingHeaderObj = ProcessingHeader()
145 145
146 146 lastBasicHeaderObj = BasicHeader(LOCALTIME)
147 147
148 148 sts = firstBasicHeaderObj.read(fp)
149 149
150 150 if not(sts):
151 151 print("[Reading] Skipping the file %s because it has not a valid header" % (filename))
152 152 return None
153 153
154 154 if not systemHeaderObj.read(fp):
155 155 return None
156 156
157 157 if not radarControllerHeaderObj.read(fp):
158 158 return None
159 159
160 160 if not processingHeaderObj.read(fp):
161 161 return None
162 162
163 163 filesize = os.path.getsize(filename)
164 164
165 165 offset = processingHeaderObj.blockSize + 24 # header size
166 166
167 167 if filesize <= offset:
168 168 print("[Reading] %s: This file has not enough data" % filename)
169 169 return None
170 170
171 171 fp.seek(-offset, 2)
172 172
173 173 sts = lastBasicHeaderObj.read(fp)
174 174
175 175 fp.close()
176 176
177 177 thisDatetime = lastBasicHeaderObj.datatime
178 178 thisTime_last_block = thisDatetime.time()
179 179
180 180 thisDatetime = firstBasicHeaderObj.datatime
181 181 thisDate = thisDatetime.date()
182 182 thisTime_first_block = thisDatetime.time()
183 183
184 184 # General case
185 185 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
186 186 #-----------o----------------------------o-----------
187 187 # startTime endTime
188 188
189 189 if endTime >= startTime:
190 190 if (thisTime_last_block < startTime) or (thisTime_first_block > endTime):
191 191 return None
192 192
193 193 return thisDatetime
194 194
195 195 # If endTime < startTime then endTime belongs to the next day
196 196
197 197 #<<<<<<<<<<<o o>>>>>>>>>>>
198 198 #-----------o----------------------------o-----------
199 199 # endTime startTime
200 200
201 201 if (thisDate == startDate) and (thisTime_last_block < startTime):
202 202 return None
203 203
204 204 if (thisDate == endDate) and (thisTime_first_block > endTime):
205 205 return None
206 206
207 207 if (thisTime_last_block < startTime) and (thisTime_first_block > endTime):
208 208 return None
209 209
210 210 return thisDatetime
211 211
212 212
213 213 def isFolderInDateRange(folder, startDate=None, endDate=None):
214 214 """
215 215 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
216 216
217 217 Inputs:
218 218 folder : nombre completo del directorio.
219 219 Su formato deberia ser "/path_root/?YYYYDDD"
220 220
221 221 siendo:
222 222 YYYY : Anio (ejemplo 2015)
223 223 DDD : Dia del anio (ejemplo 305)
224 224
225 225 startDate : fecha inicial del rango seleccionado en formato datetime.date
226 226
227 227 endDate : fecha final del rango seleccionado en formato datetime.date
228 228
229 229 Return:
230 230 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
231 231 fecha especificado, de lo contrario retorna False.
232 232 Excepciones:
233 233 Si el directorio no tiene el formato adecuado
234 234 """
235 235
236 236 basename = os.path.basename(folder)
237 237
238 238 if not isRadarFolder(basename):
239 239 print("The folder %s has not the rigth format" % folder)
240 240 return 0
241 241
242 242 if startDate and endDate:
243 243 thisDate = getDateFromRadarFolder(basename)
244 244
245 245 if thisDate < startDate:
246 246 return 0
247 247
248 248 if thisDate > endDate:
249 249 return 0
250 250
251 251 return 1
252 252
253 253
254 254 def isFileInDateRange(filename, startDate=None, endDate=None):
255 255 """
256 256 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
257 257
258 258 Inputs:
259 259 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
260 260
261 261 Su formato deberia ser "?YYYYDDDsss"
262 262
263 263 siendo:
264 264 YYYY : Anio (ejemplo 2015)
265 265 DDD : Dia del anio (ejemplo 305)
266 266 sss : set
267 267
268 268 startDate : fecha inicial del rango seleccionado en formato datetime.date
269 269
270 270 endDate : fecha final del rango seleccionado en formato datetime.date
271 271
272 272 Return:
273 273 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
274 274 fecha especificado, de lo contrario retorna False.
275 275 Excepciones:
276 276 Si el archivo no tiene el formato adecuado
277 277 """
278 278
279 279 basename = os.path.basename(filename)
280 280
281 281 if not isRadarFile(basename):
282 282 print("The filename %s has not the rigth format" % filename)
283 283 return 0
284 284
285 285 if startDate and endDate:
286 286 thisDate = getDateFromRadarFile(basename)
287 287
288 288 if thisDate < startDate:
289 289 return 0
290 290
291 291 if thisDate > endDate:
292 292 return 0
293 293
294 294 return 1
295 295
296 296
297 297 def getFileFromSet(path, ext, set):
298 298 validFilelist = []
299 299 fileList = os.listdir(path)
300 300
301 301 # 0 1234 567 89A BCDE
302 302 # H YYYY DDD SSS .ext
303 303
304 304 for thisFile in fileList:
305 305 try:
306 306 year = int(thisFile[1:5])
307 307 doy = int(thisFile[5:8])
308 308 except:
309 309 continue
310 310
311 311 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
312 312 continue
313 313
314 314 validFilelist.append(thisFile)
315 315
316 316 myfile = fnmatch.filter(
317 317 validFilelist, '*%4.4d%3.3d%3.3d*' % (year, doy, set))
318 318
319 319 if len(myfile) != 0:
320 320 return myfile[0]
321 321 else:
322 322 filename = '*%4.4d%3.3d%3.3d%s' % (year, doy, set, ext.lower())
323 323 print('the filename %s does not exist' % filename)
324 324 print('...going to the last file: ')
325 325
326 326 if validFilelist:
327 327 validFilelist = sorted(validFilelist, key=str.lower)
328 328 return validFilelist[-1]
329 329
330 330 return None
331 331
332 332
333 333 def getlastFileFromPath(path, ext):
334 334 """
335 335 Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext"
336 336 al final de la depuracion devuelve el ultimo file de la lista que quedo.
337 337
338 338 Input:
339 339 fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta
340 340 ext : extension de los files contenidos en una carpeta
341 341
342 342 Return:
343 343 El ultimo file de una determinada carpeta, no se considera el path.
344 344 """
345 345 validFilelist = []
346 346 fileList = os.listdir(path)
347 347
348 348 # 0 1234 567 89A BCDE
349 349 # H YYYY DDD SSS .ext
350 350
351 351 for thisFile in fileList:
352 352
353 353 year = thisFile[1:5]
354 354 if not isNumber(year):
355 355 continue
356 356
357 357 doy = thisFile[5:8]
358 358 if not isNumber(doy):
359 359 continue
360 360
361 361 year = int(year)
362 362 doy = int(doy)
363 363
364 364 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
365 365 continue
366 366
367 367 validFilelist.append(thisFile)
368 368
369 369 if validFilelist:
370 370 validFilelist = sorted(validFilelist, key=str.lower)
371 371 return validFilelist[-1]
372 372
373 373 return None
374 374
375 375
376 376 def isRadarFolder(folder):
377 377 try:
378 378 year = int(folder[1:5])
379 379 doy = int(folder[5:8])
380 380 except:
381 381 return 0
382 382
383 383 return 1
384 384
385 385
386 386 def isRadarFile(file):
387 387 try:
388 388 year = int(file[1:5])
389 389 doy = int(file[5:8])
390 390 set = int(file[8:11])
391 391 except:
392 392 return 0
393 393
394 394 return 1
395 395
396 396
397 397 def getDateFromRadarFile(file):
398 398 try:
399 399 year = int(file[1:5])
400 400 doy = int(file[5:8])
401 401 set = int(file[8:11])
402 402 except:
403 403 return None
404 404
405 405 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
406 406 return thisDate
407 407
408 408
409 409 def getDateFromRadarFolder(folder):
410 410 try:
411 411 year = int(folder[1:5])
412 412 doy = int(folder[5:8])
413 413 except:
414 414 return None
415 415
416 416 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
417 417 return thisDate
418 418
419 419 def parse_format(s, fmt):
420 420
421 421 for i in range(fmt.count('%')):
422 422 x = fmt.index('%')
423 423 d = DT_DIRECTIVES[fmt[x:x+2]]
424 424 fmt = fmt.replace(fmt[x:x+2], s[x:x+d])
425 425 return fmt
426 426
427 427 class Reader(object):
428 428
429 429 c = 3E8
430 430 isConfig = False
431 431 dtype = None
432 432 pathList = []
433 433 filenameList = []
434 434 datetimeList = []
435 435 filename = None
436 436 ext = None
437 437 flagIsNewFile = 1
438 438 flagDiscontinuousBlock = 0
439 439 flagIsNewBlock = 0
440 440 flagNoMoreFiles = 0
441 441 fp = None
442 442 firstHeaderSize = 0
443 443 basicHeaderSize = 24
444 444 versionFile = 1103
445 445 fileSize = None
446 446 fileSizeByHeader = None
447 447 fileIndex = -1
448 448 profileIndex = None
449 449 blockIndex = 0
450 450 nTotalBlocks = 0
451 451 maxTimeStep = 30
452 452 lastUTTime = None
453 453 datablock = None
454 454 dataOut = None
455 455 getByBlock = False
456 456 path = None
457 457 startDate = None
458 458 endDate = None
459 459 startTime = datetime.time(0, 0, 0)
460 460 endTime = datetime.time(23, 59, 59)
461 461 set = None
462 462 expLabel = ""
463 463 online = False
464 464 delay = 60
465 465 nTries = 3 # quantity tries
466 466 nFiles = 3 # number of files for searching
467 467 walk = True
468 468 getblock = False
469 469 nTxs = 1
470 470 realtime = False
471 471 blocksize = 0
472 472 blocktime = None
473 473 warnings = True
474 474 verbose = True
475 475 server = None
476 476 format = None
477 477 oneDDict = None
478 478 twoDDict = None
479 479 independentParam = None
480 480 filefmt = None
481 481 folderfmt = None
482 482 open_file = open
483 483 open_mode = 'rb'
484 484
485 485 def run(self):
486 486
487 487 raise NotImplementedError
488 488
489 489 def getAllowedArgs(self):
490 490 if hasattr(self, '__attrs__'):
491 491 return self.__attrs__
492 492 else:
493 493 return inspect.getargspec(self.run).args
494 494
495 495 def set_kwargs(self, **kwargs):
496 496
497 497 for key, value in kwargs.items():
498 498 setattr(self, key, value)
499 499
500 500 def find_folders(self, path, startDate, endDate, folderfmt, last=False):
501 501
502 502 folders = [x for f in path.split(',')
503 503 for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
504 504 folders.sort()
505 505
506 506 if last:
507 507 folders = [folders[-1]]
508 508
509 509 for folder in folders:
510 510 try:
511 511 dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
512 512 if dt >= startDate and dt <= endDate:
513 513 yield os.path.join(path, folder)
514 514 else:
515 515 log.log('Skiping folder {}'.format(folder), self.name)
516 516 except Exception as e:
517 517 log.log('Skiping folder {}'.format(folder), self.name)
518 518 continue
519 519 return
520 520
521 521 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
522 522 expLabel='', last=False):
523 523
524 524 for path in folders:
525 525 files = glob.glob1(path, '*{}'.format(ext))
526 526 files.sort()
527 527 if last:
528 528 if files:
529 529 fo = files[-1]
530 530 try:
531 531 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
532 532 yield os.path.join(path, expLabel, fo)
533 533 except Exception as e:
534 534 pass
535 535 return
536 536 else:
537 537 return
538 538
539 539 for fo in files:
540 540 try:
541 541 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
542 542 if dt >= startDate and dt <= endDate:
543 543 yield os.path.join(path, expLabel, fo)
544 544 else:
545 545 log.log('Skiping file {}'.format(fo), self.name)
546 546 except Exception as e:
547 547 log.log('Skiping file {}'.format(fo), self.name)
548 548 continue
549 549
550 550 def searchFilesOffLine(self, path, startDate, endDate,
551 551 expLabel, ext, walk,
552 552 filefmt, folderfmt):
553 553 """Search files in offline mode for the given arguments
554 554
555 555 Return:
556 556 Generator of files
557 557 """
558 558
559 559 if walk:
560 560 folders = self.find_folders(
561 561 path, startDate, endDate, folderfmt)
562 562 else:
563 563 folders = path.split(',')
564 564
565 565 return self.find_files(
566 566 folders, ext, filefmt, startDate, endDate, expLabel)
567 567
568 568 def searchFilesOnLine(self, path, startDate, endDate,
569 569 expLabel, ext, walk,
570 570 filefmt, folderfmt):
571 571 """Search for the last file of the last folder
572 572
573 573 Arguments:
574 574 path : carpeta donde estan contenidos los files que contiene data
575 575 expLabel : Nombre del subexperimento (subfolder)
576 576 ext : extension de los files
577 577 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
578 578
579 579 Return:
580 580 generator with the full path of last filename
581 581 """
582 582
583 583 if walk:
584 584 folders = self.find_folders(
585 585 path, startDate, endDate, folderfmt, last=True)
586 586 else:
587 587 folders = path.split(',')
588 588
589 589 return self.find_files(
590 590 folders, ext, filefmt, startDate, endDate, expLabel, last=True)
591 591
592 592 def setNextFile(self):
593 593 """Set the next file to be readed open it and parse de file header"""
594 594
595 595 while True:
596 596 if self.fp != None:
597 597 self.fp.close()
598 598
599 599 if self.online:
600 600 newFile = self.setNextFileOnline()
601 601 else:
602 602 newFile = self.setNextFileOffline()
603 603
604 604 if not(newFile):
605 605 if self.online:
606 606 raise schainpy.admin.SchainError('Time to wait for new files reach')
607 607 else:
608 608 if self.fileIndex == -1:
609 609 raise schainpy.admin.SchainWarning('No files found in the given path')
610 610 else:
611 611 raise schainpy.admin.SchainWarning('No more files to read')
612 612
613 613 if self.verifyFile(self.filename):
614 614 break
615 615
616 616 log.log('Opening file: %s' % self.filename, self.name)
617 617
618 618 self.readFirstHeader()
619 619 self.nReadBlocks = 0
620 620
621 621 def setNextFileOnline(self):
622 622 """Check for the next file to be readed in online mode.
623 623
624 624 Set:
625 625 self.filename
626 626 self.fp
627 627 self.filesize
628 628
629 629 Return:
630 630 boolean
631 631
632 632 """
633 633 nextFile = True
634 634 nextDay = False
635 635
636 636 for nFiles in range(self.nFiles+1):
637 637 for nTries in range(self.nTries):
638 638 fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
639 639 if fullfilename is not None:
640 640 break
641 641 log.warning(
642 642 "Waiting %0.2f sec for the next file: \"%s\" , try %02d ..." % (self.delay, filename, nTries + 1),
643 643 self.name)
644 644 time.sleep(self.delay)
645 645 nextFile = False
646 646 continue
647 647
648 648 if fullfilename is not None:
649 649 break
650 650
651 651 self.nTries = 1
652 652 nextFile = True
653 653
654 654 if nFiles == (self.nFiles - 1):
655 655 log.log('Trying with next day...', self.name)
656 656 nextDay = True
657 657 self.nTries = 3
658 658
659 659 if fullfilename:
660 660 self.fileSize = os.path.getsize(fullfilename)
661 661 self.filename = fullfilename
662 662 self.flagIsNewFile = 1
663 663 if self.fp != None:
664 664 self.fp.close()
665 665 self.fp = self.open_file(fullfilename, self.open_mode)
666 666 self.flagNoMoreFiles = 0
667 667 self.fileIndex += 1
668 668 return 1
669 669 else:
670 670 return 0
671 671
672 672 def setNextFileOffline(self):
673 673 """Open the next file to be readed in offline mode"""
674 674
675 675 try:
676 676 filename = next(self.filenameList)
677 677 self.fileIndex +=1
678 678 except StopIteration:
679 679 self.flagNoMoreFiles = 1
680 680 return 0
681 681
682 682 self.filename = filename
683 683 self.fileSize = os.path.getsize(filename)
684 684 self.fp = self.open_file(filename, self.open_mode)
685 685 self.flagIsNewFile = 1
686 686
687 687 return 1
688 688
689 689 @staticmethod
690 690 def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
691 691 """Check if the given datetime is in range"""
692
693 if startDate <= dt.date() <= endDate:
694 if startTime <= dt.time() <= endTime:
692 startDateTime= datetime.datetime.combine(startDate,startTime)
693 endDateTime = datetime.datetime.combine(endDate,endTime)
694 if startDateTime <= dt <= endDateTime:
695 695 return True
696 696 return False
697 697
698 698 def verifyFile(self, filename):
699 699 """Check for a valid file
700 700
701 701 Arguments:
702 702 filename -- full path filename
703 703
704 704 Return:
705 705 boolean
706 706 """
707 707
708 708 return True
709 709
710 710 def checkForRealPath(self, nextFile, nextDay):
711 711 """Check if the next file to be readed exists"""
712 712
713 713 raise NotImplementedError
714 714
715 715 def readFirstHeader(self):
716 716 """Parse the file header"""
717 717
718 718 pass
719 719
720 720 def waitDataBlock(self, pointer_location, blocksize=None):
721 721 """
722 722 """
723 723
724 724 currentPointer = pointer_location
725 725 if blocksize is None:
726 726 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
727 727 else:
728 728 neededSize = blocksize
729 729
730 730 for nTries in range(self.nTries):
731 731 self.fp.close()
732 732 self.fp = open(self.filename, 'rb')
733 733 self.fp.seek(currentPointer)
734 734
735 735 self.fileSize = os.path.getsize(self.filename)
736 736 currentSize = self.fileSize - currentPointer
737 737
738 738 if (currentSize >= neededSize):
739 739 return 1
740 740
741 741 log.warning(
742 742 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
743 743 self.name
744 744 )
745 745 time.sleep(self.delay)
746 746
747 747 return 0
748 748
749 749 class JRODataReader(Reader):
750 750
751 751 utc = 0
752 752 nReadBlocks = 0
753 753 foldercounter = 0
754 754 firstHeaderSize = 0
755 755 basicHeaderSize = 24
756 756 __isFirstTimeOnline = 1
757 757 filefmt = "*%Y%j***"
758 758 folderfmt = "*%Y%j"
759 759 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'online', 'delay', 'walk']
760 760
761 761 def getDtypeWidth(self):
762 762
763 763 dtype_index = get_dtype_index(self.dtype)
764 764 dtype_width = get_dtype_width(dtype_index)
765 765
766 766 return dtype_width
767 767
768 768 def checkForRealPath(self, nextFile, nextDay):
769 769 """Check if the next file to be readed exists.
770 770
771 771 Example :
772 772 nombre correcto del file es .../.../D2009307/P2009307367.ext
773 773
774 774 Entonces la funcion prueba con las siguientes combinaciones
775 775 .../.../y2009307367.ext
776 776 .../.../Y2009307367.ext
777 777 .../.../x2009307/y2009307367.ext
778 778 .../.../x2009307/Y2009307367.ext
779 779 .../.../X2009307/y2009307367.ext
780 780 .../.../X2009307/Y2009307367.ext
781 781 siendo para este caso, la ultima combinacion de letras, identica al file buscado
782 782
783 783 Return:
784 784 str -- fullpath of the file
785 785 """
786 786
787 787
788 788 if nextFile:
789 789 self.set += 1
790 790 if nextDay:
791 791 self.set = 0
792 792 self.doy += 1
793 793 foldercounter = 0
794 794 prefixDirList = [None, 'd', 'D']
795 795 if self.ext.lower() == ".r": # voltage
796 796 prefixFileList = ['d', 'D']
797 797 elif self.ext.lower() == ".pdata": # spectra
798 798 prefixFileList = ['p', 'P']
799 799
800 800 # barrido por las combinaciones posibles
801 801 for prefixDir in prefixDirList:
802 802 thispath = self.path
803 803 if prefixDir != None:
804 804 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
805 805 if foldercounter == 0:
806 806 thispath = os.path.join(self.path, "%s%04d%03d" %
807 807 (prefixDir, self.year, self.doy))
808 808 else:
809 809 thispath = os.path.join(self.path, "%s%04d%03d_%02d" % (
810 810 prefixDir, self.year, self.doy, foldercounter))
811 811 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
812 812 # formo el nombre del file xYYYYDDDSSS.ext
813 813 filename = "%s%04d%03d%03d%s" % (prefixFile, self.year, self.doy, self.set, self.ext)
814 814 fullfilename = os.path.join(
815 815 thispath, filename)
816 816
817 817 if os.path.exists(fullfilename):
818 818 return fullfilename, filename
819 819
820 820 return None, filename
821 821
822 822 def __waitNewBlock(self):
823 823 """
824 824 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
825 825
826 826 Si el modo de lectura es OffLine siempre retorn 0
827 827 """
828 828 if not self.online:
829 829 return 0
830 830
831 831 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
832 832 return 0
833 833
834 834 currentPointer = self.fp.tell()
835 835
836 836 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
837 837
838 838 for nTries in range(self.nTries):
839 839
840 840 self.fp.close()
841 841 self.fp = open(self.filename, 'rb')
842 842 self.fp.seek(currentPointer)
843 843
844 844 self.fileSize = os.path.getsize(self.filename)
845 845 currentSize = self.fileSize - currentPointer
846 846
847 847 if (currentSize >= neededSize):
848 848 self.basicHeaderObj.read(self.fp)
849 849 return 1
850 850
851 851 if self.fileSize == self.fileSizeByHeader:
852 852 # self.flagEoF = True
853 853 return 0
854 854
855 855 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
856 856 time.sleep(self.delay)
857 857
858 858 return 0
859 859
860 860 def __setNewBlock(self):
861 861
862 862 if self.fp == None:
863 863 return 0
864 864
865 865 if self.flagIsNewFile:
866 866 self.lastUTTime = self.basicHeaderObj.utc
867 867 return 1
868 868
869 869 if self.realtime:
870 870 self.flagDiscontinuousBlock = 1
871 871 if not(self.setNextFile()):
872 872 return 0
873 873 else:
874 874 return 1
875 875
876 876 currentSize = self.fileSize - self.fp.tell()
877 877 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
878 878
879 879 if (currentSize >= neededSize):
880 880 self.basicHeaderObj.read(self.fp)
881 881 self.lastUTTime = self.basicHeaderObj.utc
882 882 return 1
883 883
884 884 if self.__waitNewBlock():
885 885 self.lastUTTime = self.basicHeaderObj.utc
886 886 return 1
887 887
888 888 if not(self.setNextFile()):
889 889 return 0
890 890
891 891 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
892 892 self.lastUTTime = self.basicHeaderObj.utc
893 893
894 894 self.flagDiscontinuousBlock = 0
895 895
896 896 if deltaTime > self.maxTimeStep:
897 897 self.flagDiscontinuousBlock = 1
898 898
899 899 return 1
900 900
901 901 def readNextBlock(self):
902 902
903 903 while True:
904 904 if not(self.__setNewBlock()):
905 905 continue
906 906
907 907 if not(self.readBlock()):
908 908 return 0
909 909
910 910 self.getBasicHeader()
911 911
912 912 if not self.isDateTimeInRange(self.dataOut.datatime, self.startDate, self.endDate, self.startTime, self.endTime):
913 913 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
914 914 self.processingHeaderObj.dataBlocksPerFile,
915 915 self.dataOut.datatime.ctime()))
916 916 continue
917 917
918 918 break
919 919
920 920 if self.verbose:
921 921 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
922 922 self.processingHeaderObj.dataBlocksPerFile,
923 923 self.dataOut.datatime.ctime()))
924 924 return 1
925 925
926 926 def readFirstHeader(self):
927 927
928 928 self.basicHeaderObj.read(self.fp)
929 929 self.systemHeaderObj.read(self.fp)
930 930 self.radarControllerHeaderObj.read(self.fp)
931 931 self.processingHeaderObj.read(self.fp)
932 932 self.firstHeaderSize = self.basicHeaderObj.size
933 933
934 934 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
935 935 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
936 936 if datatype == 0:
937 937 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
938 938 elif datatype == 1:
939 939 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
940 940 elif datatype == 2:
941 941 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
942 942 elif datatype == 3:
943 943 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
944 944 elif datatype == 4:
945 945 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
946 946 elif datatype == 5:
947 947 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
948 948 else:
949 949 raise ValueError('Data type was not defined')
950 950
951 951 self.dtype = datatype_str
952 952 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
953 953 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
954 954 self.firstHeaderSize + self.basicHeaderSize * \
955 955 (self.processingHeaderObj.dataBlocksPerFile - 1)
956 956 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
957 957 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
958 958 self.getBlockDimension()
959 959
960 960 def verifyFile(self, filename):
961 961
962 962 flag = True
963 963
964 964 try:
965 965 fp = open(filename, 'rb')
966 966 except IOError:
967 967 log.error("File {} can't be opened".format(filename), self.name)
968 968 return False
969 969
970 970 if self.online and self.waitDataBlock(0):
971 971 pass
972 972
973 973 basicHeaderObj = BasicHeader(LOCALTIME)
974 974 systemHeaderObj = SystemHeader()
975 975 radarControllerHeaderObj = RadarControllerHeader()
976 976 processingHeaderObj = ProcessingHeader()
977 977
978 978 if not(basicHeaderObj.read(fp)):
979 979 flag = False
980 980 if not(systemHeaderObj.read(fp)):
981 981 flag = False
982 982 if not(radarControllerHeaderObj.read(fp)):
983 983 flag = False
984 984 if not(processingHeaderObj.read(fp)):
985 985 flag = False
986 986 if not self.online:
987 987 dt1 = basicHeaderObj.datatime
988 988 pos = self.fileSize-processingHeaderObj.blockSize-24
989 989 if pos<0:
990 990 flag = False
991 991 log.error('Invalid size for file: {}'.format(self.filename), self.name)
992 992 else:
993 993 fp.seek(pos)
994 994 if not(basicHeaderObj.read(fp)):
995 995 flag = False
996 996 dt2 = basicHeaderObj.datatime
997 997 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
998 998 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
999 999 flag = False
1000 1000
1001 1001 fp.close()
1002 1002 return flag
1003 1003
1004 1004 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1005 1005
1006 1006 path_empty = True
1007 1007
1008 1008 dateList = []
1009 1009 pathList = []
1010 1010
1011 1011 multi_path = path.split(',')
1012 1012
1013 1013 if not walk:
1014 1014
1015 1015 for single_path in multi_path:
1016 1016
1017 1017 if not os.path.isdir(single_path):
1018 1018 continue
1019 1019
1020 1020 fileList = glob.glob1(single_path, "*" + ext)
1021 1021
1022 1022 if not fileList:
1023 1023 continue
1024 1024
1025 1025 path_empty = False
1026 1026
1027 1027 fileList.sort()
1028 1028
1029 1029 for thisFile in fileList:
1030 1030
1031 1031 if not os.path.isfile(os.path.join(single_path, thisFile)):
1032 1032 continue
1033 1033
1034 1034 if not isRadarFile(thisFile):
1035 1035 continue
1036 1036
1037 1037 if not isFileInDateRange(thisFile, startDate, endDate):
1038 1038 continue
1039 1039
1040 1040 thisDate = getDateFromRadarFile(thisFile)
1041 1041
1042 1042 if thisDate in dateList or single_path in pathList:
1043 1043 continue
1044 1044
1045 1045 dateList.append(thisDate)
1046 1046 pathList.append(single_path)
1047 1047
1048 1048 else:
1049 1049 for single_path in multi_path:
1050 1050
1051 1051 if not os.path.isdir(single_path):
1052 1052 continue
1053 1053
1054 1054 dirList = []
1055 1055
1056 1056 for thisPath in os.listdir(single_path):
1057 1057
1058 1058 if not os.path.isdir(os.path.join(single_path, thisPath)):
1059 1059 continue
1060 1060
1061 1061 if not isRadarFolder(thisPath):
1062 1062 continue
1063 1063
1064 1064 if not isFolderInDateRange(thisPath, startDate, endDate):
1065 1065 continue
1066 1066
1067 1067 dirList.append(thisPath)
1068 1068
1069 1069 if not dirList:
1070 1070 continue
1071 1071
1072 1072 dirList.sort()
1073 1073
1074 1074 for thisDir in dirList:
1075 1075
1076 1076 datapath = os.path.join(single_path, thisDir, expLabel)
1077 1077 fileList = glob.glob1(datapath, "*" + ext)
1078 1078
1079 1079 if not fileList:
1080 1080 continue
1081 1081
1082 1082 path_empty = False
1083 1083
1084 1084 thisDate = getDateFromRadarFolder(thisDir)
1085 1085
1086 1086 pathList.append(datapath)
1087 1087 dateList.append(thisDate)
1088 1088
1089 1089 dateList.sort()
1090 1090
1091 1091 if walk:
1092 1092 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1093 1093 else:
1094 1094 pattern_path = multi_path[0]
1095 1095
1096 1096 if path_empty:
1097 1097 raise schainpy.admin.SchainError("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1098 1098 else:
1099 1099 if not dateList:
1100 1100 raise schainpy.admin.SchainError("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1101 1101
1102 1102 if include_path:
1103 1103 return dateList, pathList
1104 1104
1105 1105 return dateList
1106 1106
1107 1107 def setup(self, **kwargs):
1108 1108
1109 1109 self.set_kwargs(**kwargs)
1110 1110 if not self.ext.startswith('.'):
1111 1111 self.ext = '.{}'.format(self.ext)
1112 1112
1113 1113 if self.server is not None:
1114 1114 if 'tcp://' in self.server:
1115 1115 address = server
1116 1116 else:
1117 1117 address = 'ipc:///tmp/%s' % self.server
1118 1118 self.server = address
1119 1119 self.context = zmq.Context()
1120 1120 self.receiver = self.context.socket(zmq.PULL)
1121 1121 self.receiver.connect(self.server)
1122 1122 time.sleep(0.5)
1123 1123 print('[Starting] ReceiverData from {}'.format(self.server))
1124 1124 else:
1125 1125 self.server = None
1126 1126 if self.path == None:
1127 1127 raise ValueError("[Reading] The path is not valid")
1128 1128
1129 1129 if self.online:
1130 1130 log.log("[Reading] Searching files in online mode...", self.name)
1131 1131
1132 1132 for nTries in range(self.nTries):
1133 1133 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1134 1134 self.endDate, self.expLabel, self.ext, self.walk,
1135 1135 self.filefmt, self.folderfmt)
1136 1136
1137 1137 try:
1138 1138 fullpath = next(fullpath)
1139 1139 except:
1140 1140 fullpath = None
1141 1141
1142 1142 if fullpath:
1143 1143 break
1144 1144
1145 1145 log.warning(
1146 1146 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1147 1147 self.delay, self.path, nTries + 1),
1148 1148 self.name)
1149 1149 time.sleep(self.delay)
1150 1150
1151 1151 if not(fullpath):
1152 1152 raise schainpy.admin.SchainError(
1153 1153 'There isn\'t any valid file in {}'.format(self.path))
1154 1154
1155 1155 pathname, filename = os.path.split(fullpath)
1156 1156 self.year = int(filename[1:5])
1157 1157 self.doy = int(filename[5:8])
1158 1158 self.set = int(filename[8:11]) - 1
1159 1159 else:
1160 1160 log.log("Searching files in {}".format(self.path), self.name)
1161 1161 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1162 1162 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1163 1163
1164 1164 self.setNextFile()
1165 1165
1166 1166 return
1167 1167
1168 1168 def getBasicHeader(self):
1169 1169
1170 1170 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1171 1171 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1172 1172
1173 1173 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1174 1174
1175 1175 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1176 1176
1177 1177 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1178 1178
1179 1179 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1180 1180
1181 1181 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1182 1182
1183 1183 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1184 1184
1185 1185 def getFirstHeader(self):
1186 1186
1187 1187 raise NotImplementedError
1188 1188
1189 1189 def getData(self):
1190 1190
1191 1191 raise NotImplementedError
1192 1192
1193 1193 def hasNotDataInBuffer(self):
1194 1194
1195 1195 raise NotImplementedError
1196 1196
1197 1197 def readBlock(self):
1198 1198
1199 1199 raise NotImplementedError
1200 1200
1201 1201 def isEndProcess(self):
1202 1202
1203 1203 return self.flagNoMoreFiles
1204 1204
1205 1205 def printReadBlocks(self):
1206 1206
1207 1207 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1208 1208
1209 1209 def printTotalBlocks(self):
1210 1210
1211 1211 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1212 1212
1213 1213 def run(self, **kwargs):
1214 1214 """
1215 1215
1216 1216 Arguments:
1217 1217 path :
1218 1218 startDate :
1219 1219 endDate :
1220 1220 startTime :
1221 1221 endTime :
1222 1222 set :
1223 1223 expLabel :
1224 1224 ext :
1225 1225 online :
1226 1226 delay :
1227 1227 walk :
1228 1228 getblock :
1229 1229 nTxs :
1230 1230 realtime :
1231 1231 blocksize :
1232 1232 blocktime :
1233 1233 skip :
1234 1234 cursor :
1235 1235 warnings :
1236 1236 server :
1237 1237 verbose :
1238 1238 format :
1239 1239 oneDDict :
1240 1240 twoDDict :
1241 1241 independentParam :
1242 1242 """
1243 1243
1244 1244 if not(self.isConfig):
1245 1245 self.setup(**kwargs)
1246 1246 self.isConfig = True
1247 1247 if self.server is None:
1248 1248 self.getData()
1249 1249 else:
1250 1250 self.getFromServer()
1251 1251
1252 1252
1253 1253 class JRODataWriter(Reader):
1254 1254
1255 1255 """
1256 1256 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1257 1257 de los datos siempre se realiza por bloques.
1258 1258 """
1259 1259
1260 1260 setFile = None
1261 1261 profilesPerBlock = None
1262 1262 blocksPerFile = None
1263 1263 nWriteBlocks = 0
1264 1264 fileDate = None
1265 1265
1266 1266 def __init__(self, dataOut=None):
1267 1267 raise NotImplementedError
1268 1268
1269 1269 def hasAllDataInBuffer(self):
1270 1270 raise NotImplementedError
1271 1271
1272 1272 def setBlockDimension(self):
1273 1273 raise NotImplementedError
1274 1274
1275 1275 def writeBlock(self):
1276 1276 raise NotImplementedError
1277 1277
1278 1278 def putData(self):
1279 1279 raise NotImplementedError
1280 1280
1281 1281 def getDtypeWidth(self):
1282 1282
1283 1283 dtype_index = get_dtype_index(self.dtype)
1284 1284 dtype_width = get_dtype_width(dtype_index)
1285 1285
1286 1286 return dtype_width
1287 1287
1288 1288 def getProcessFlags(self):
1289 1289
1290 1290 processFlags = 0
1291 1291
1292 1292 dtype_index = get_dtype_index(self.dtype)
1293 1293 procflag_dtype = get_procflag_dtype(dtype_index)
1294 1294
1295 1295 processFlags += procflag_dtype
1296 1296
1297 1297 if self.dataOut.flagDecodeData:
1298 1298 processFlags += PROCFLAG.DECODE_DATA
1299 1299
1300 1300 if self.dataOut.flagDeflipData:
1301 1301 processFlags += PROCFLAG.DEFLIP_DATA
1302 1302
1303 1303 if self.dataOut.code is not None:
1304 1304 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1305 1305
1306 1306 if self.dataOut.nCohInt > 1:
1307 1307 processFlags += PROCFLAG.COHERENT_INTEGRATION
1308 1308
1309 1309 if self.dataOut.type == "Spectra":
1310 1310 if self.dataOut.nIncohInt > 1:
1311 1311 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1312 1312
1313 1313 if self.dataOut.data_dc is not None:
1314 1314 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1315 1315
1316 1316 if self.dataOut.flagShiftFFT:
1317 1317 processFlags += PROCFLAG.SHIFT_FFT_DATA
1318 1318
1319 1319 return processFlags
1320 1320
1321 1321 def setBasicHeader(self):
1322 1322
1323 1323 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1324 1324 self.basicHeaderObj.version = self.versionFile
1325 1325 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1326 1326 utc = numpy.floor(self.dataOut.utctime)
1327 1327 milisecond = (self.dataOut.utctime - utc) * 1000.0
1328 1328 self.basicHeaderObj.utc = utc
1329 1329 self.basicHeaderObj.miliSecond = milisecond
1330 1330 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1331 1331 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1332 1332 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1333 1333
1334 1334 def setFirstHeader(self):
1335 1335 """
1336 1336 Obtiene una copia del First Header
1337 1337
1338 1338 Affected:
1339 1339
1340 1340 self.basicHeaderObj
1341 1341 self.systemHeaderObj
1342 1342 self.radarControllerHeaderObj
1343 1343 self.processingHeaderObj self.
1344 1344
1345 1345 Return:
1346 1346 None
1347 1347 """
1348 1348
1349 1349 raise NotImplementedError
1350 1350
1351 1351 def __writeFirstHeader(self):
1352 1352 """
1353 1353 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1354 1354
1355 1355 Affected:
1356 1356 __dataType
1357 1357
1358 1358 Return:
1359 1359 None
1360 1360 """
1361 1361
1362 1362 # CALCULAR PARAMETROS
1363 1363
1364 1364 sizeLongHeader = self.systemHeaderObj.size + \
1365 1365 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1366 1366 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1367 1367
1368 1368 self.basicHeaderObj.write(self.fp)
1369 1369 self.systemHeaderObj.write(self.fp)
1370 1370 self.radarControllerHeaderObj.write(self.fp)
1371 1371 self.processingHeaderObj.write(self.fp)
1372 1372
1373 1373 def __setNewBlock(self):
1374 1374 """
1375 1375 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1376 1376
1377 1377 Return:
1378 1378 0 : si no pudo escribir nada
1379 1379 1 : Si escribio el Basic el First Header
1380 1380 """
1381 1381 if self.fp == None:
1382 1382 self.setNextFile()
1383 1383
1384 1384 if self.flagIsNewFile:
1385 1385 return 1
1386 1386
1387 1387 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1388 1388 self.basicHeaderObj.write(self.fp)
1389 1389 return 1
1390 1390
1391 1391 if not(self.setNextFile()):
1392 1392 return 0
1393 1393
1394 1394 return 1
1395 1395
1396 1396 def writeNextBlock(self):
1397 1397 """
1398 1398 Selecciona el bloque siguiente de datos y los escribe en un file
1399 1399
1400 1400 Return:
1401 1401 0 : Si no hizo pudo escribir el bloque de datos
1402 1402 1 : Si no pudo escribir el bloque de datos
1403 1403 """
1404 1404 if not(self.__setNewBlock()):
1405 1405 return 0
1406 1406
1407 1407 self.writeBlock()
1408 1408
1409 1409 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1410 1410 self.processingHeaderObj.dataBlocksPerFile))
1411 1411
1412 1412 return 1
1413 1413
1414 1414 def setNextFile(self):
1415 1415 """Determina el siguiente file que sera escrito
1416 1416
1417 1417 Affected:
1418 1418 self.filename
1419 1419 self.subfolder
1420 1420 self.fp
1421 1421 self.setFile
1422 1422 self.flagIsNewFile
1423 1423
1424 1424 Return:
1425 1425 0 : Si el archivo no puede ser escrito
1426 1426 1 : Si el archivo esta listo para ser escrito
1427 1427 """
1428 1428 ext = self.ext
1429 1429 path = self.path
1430 1430
1431 1431 if self.fp != None:
1432 1432 self.fp.close()
1433 1433
1434 1434 if not os.path.exists(path):
1435 1435 os.mkdir(path)
1436 1436
1437 1437 timeTuple = time.localtime(self.dataOut.utctime)
1438 1438 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1439 1439
1440 1440 fullpath = os.path.join(path, subfolder)
1441 1441 setFile = self.setFile
1442 1442
1443 1443 if not(os.path.exists(fullpath)):
1444 1444 os.mkdir(fullpath)
1445 1445 setFile = -1 # inicializo mi contador de seteo
1446 1446 else:
1447 1447 filesList = os.listdir(fullpath)
1448 1448 if len(filesList) > 0:
1449 1449 filesList = sorted(filesList, key=str.lower)
1450 1450 filen = filesList[-1]
1451 1451 # el filename debera tener el siguiente formato
1452 1452 # 0 1234 567 89A BCDE (hex)
1453 1453 # x YYYY DDD SSS .ext
1454 1454 if isNumber(filen[8:11]):
1455 1455 # inicializo mi contador de seteo al seteo del ultimo file
1456 1456 setFile = int(filen[8:11])
1457 1457 else:
1458 1458 setFile = -1
1459 1459 else:
1460 1460 setFile = -1 # inicializo mi contador de seteo
1461 1461
1462 1462 setFile += 1
1463 1463
1464 1464 # If this is a new day it resets some values
1465 1465 if self.dataOut.datatime.date() > self.fileDate:
1466 1466 setFile = 0
1467 1467 self.nTotalBlocks = 0
1468 1468
1469 1469 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1470 1470 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1471 1471
1472 1472 filename = os.path.join(path, subfolder, filen)
1473 1473
1474 1474 fp = open(filename, 'wb')
1475 1475
1476 1476 self.blockIndex = 0
1477 1477 self.filename = filename
1478 1478 self.subfolder = subfolder
1479 1479 self.fp = fp
1480 1480 self.setFile = setFile
1481 1481 self.flagIsNewFile = 1
1482 1482 self.fileDate = self.dataOut.datatime.date()
1483 1483 self.setFirstHeader()
1484 1484
1485 1485 print('[Writing] Opening file: %s' % self.filename)
1486 1486
1487 1487 self.__writeFirstHeader()
1488 1488
1489 1489 return 1
1490 1490
1491 1491 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1492 1492 """
1493 1493 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1494 1494
1495 1495 Inputs:
1496 1496 path : directory where data will be saved
1497 1497 profilesPerBlock : number of profiles per block
1498 1498 set : initial file set
1499 1499 datatype : An integer number that defines data type:
1500 1500 0 : int8 (1 byte)
1501 1501 1 : int16 (2 bytes)
1502 1502 2 : int32 (4 bytes)
1503 1503 3 : int64 (8 bytes)
1504 1504 4 : float32 (4 bytes)
1505 1505 5 : double64 (8 bytes)
1506 1506
1507 1507 Return:
1508 1508 0 : Si no realizo un buen seteo
1509 1509 1 : Si realizo un buen seteo
1510 1510 """
1511 1511
1512 1512 if ext == None:
1513 1513 ext = self.ext
1514 1514
1515 1515 self.ext = ext.lower()
1516 1516
1517 1517 self.path = path
1518 1518
1519 1519 if set is None:
1520 1520 self.setFile = -1
1521 1521 else:
1522 1522 self.setFile = set - 1
1523 1523
1524 1524 self.blocksPerFile = blocksPerFile
1525 1525 self.profilesPerBlock = profilesPerBlock
1526 1526 self.dataOut = dataOut
1527 1527 self.fileDate = self.dataOut.datatime.date()
1528 1528 self.dtype = self.dataOut.dtype
1529 1529
1530 1530 if datatype is not None:
1531 1531 self.dtype = get_numpy_dtype(datatype)
1532 1532
1533 1533 if not(self.setNextFile()):
1534 1534 print("[Writing] There isn't a next file")
1535 1535 return 0
1536 1536
1537 1537 self.setBlockDimension()
1538 1538
1539 1539 return 1
1540 1540
1541 1541 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1542 1542
1543 1543 if not(self.isConfig):
1544 1544
1545 1545 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1546 1546 set=set, ext=ext, datatype=datatype, **kwargs)
1547 1547 self.isConfig = True
1548 1548
1549 1549 self.dataOut = dataOut
1550 1550 self.putData()
1551 1551 return self.dataOut
1552 1552
1553 1553 @MPDecorator
1554 1554 class printInfo(Operation):
1555 1555
1556 1556 def __init__(self):
1557 1557
1558 1558 Operation.__init__(self)
1559 1559 self.__printInfo = True
1560 1560
1561 1561 def run(self, dataOut, headers = ['systemHeaderObj', 'radarControllerHeaderObj', 'processingHeaderObj']):
1562 1562 if self.__printInfo == False:
1563 1563 return
1564 1564
1565 1565 for header in headers:
1566 1566 if hasattr(dataOut, header):
1567 1567 obj = getattr(dataOut, header)
1568 1568 if hasattr(obj, 'printInfo'):
1569 1569 obj.printInfo()
1570 1570 else:
1571 1571 print(obj)
1572 1572 else:
1573 1573 log.warning('Header {} Not found in object'.format(header))
1574 1574
1575 1575 self.__printInfo = False
Show file before | Show file after | Hide comments
@@ -1,662 +1,664
1 1 '''
2 2 Created on Set 9, 2015
3 3
4 4 @author: roj-idl71 Karim Kuyeng
5 5
6 6 @update: 2021, Joab Apaza
7 7 '''
8 8
9 9 import os
10 10 import sys
11 11 import glob
12 12 import fnmatch
13 13 import datetime
14 14 import time
15 15 import re
16 16 import h5py
17 17 import numpy
18 18
19 19 try:
20 20 from gevent import sleep
21 21 except:
22 22 from time import sleep
23 23
24 24 from schainpy.model.data.jroheaderIO import RadarControllerHeader, SystemHeader
25 25 from schainpy.model.data.jrodata import Voltage
26 26 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
27 27 from numpy import imag
28 28
29 29
30 30 class AMISRReader(ProcessingUnit):
31 31 '''
32 32 classdocs
33 33 '''
34 34
35 35 def __init__(self):
36 36 '''
37 37 Constructor
38 38 '''
39 39
40 40 ProcessingUnit.__init__(self)
41 41
42 42 self.set = None
43 43 self.subset = None
44 44 self.extension_file = '.h5'
45 45 self.dtc_str = 'dtc'
46 46 self.dtc_id = 0
47 47 self.status = True
48 48 self.isConfig = False
49 49 self.dirnameList = []
50 50 self.filenameList = []
51 51 self.fileIndex = None
52 52 self.flagNoMoreFiles = False
53 53 self.flagIsNewFile = 0
54 54 self.filename = ''
55 55 self.amisrFilePointer = None
56 56 self.realBeamCode = []
57 57
58 58 self.dataShape = None
59 59
60 60
61 61
62 62 self.profileIndex = 0
63 63
64 64
65 65 self.beamCodeByFrame = None
66 66 self.radacTimeByFrame = None
67 67
68 68 self.dataset = None
69 69
70 70
71 71
72 72
73 73 self.__firstFile = True
74 74
75 75 self.buffer = None
76 76
77 77
78 78 self.timezone = 'ut'
79 79
80 80 self.__waitForNewFile = 20
81 81 self.__filename_online = None
82 82 #Is really necessary create the output object in the initializer
83 83 self.dataOut = Voltage()
84 84 self.dataOut.error=False
85 85
86 86
87 87 def setup(self,path=None,
88 88 startDate=None,
89 89 endDate=None,
90 90 startTime=None,
91 91 endTime=None,
92 92 walk=True,
93 93 timezone='ut',
94 94 all=0,
95 95 code = None,
96 96 nCode = 0,
97 97 nBaud = 0,
98 98 online=False):
99 99
100 100
101 101
102 102 self.timezone = timezone
103 103 self.all = all
104 104 self.online = online
105 105
106 106 self.code = code
107 107 self.nCode = int(nCode)
108 108 self.nBaud = int(nBaud)
109 109
110 110
111 111
112 112 #self.findFiles()
113 113 if not(online):
114 114 #Busqueda de archivos offline
115 115 self.searchFilesOffLine(path, startDate, endDate, startTime, endTime, walk)
116 116 else:
117 117 self.searchFilesOnLine(path, startDate, endDate, startTime,endTime,walk)
118 118
119 119 if not(self.filenameList):
120 120 print("There is no files into the folder: %s"%(path))
121 121 sys.exit()
122 122
123 123 self.fileIndex = 0
124 124
125 125 self.readNextFile(online)
126 126
127 127 '''
128 128 Add code
129 129 '''
130 130 self.isConfig = True
131 131 # print("Setup Done")
132 132 pass
133 133
134 134
135 135 def readAMISRHeader(self,fp):
136 136
137 137 if self.isConfig and (not self.flagNoMoreFiles):
138 138 newShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
139 139 if self.dataShape != newShape and newShape != None:
140 140 print("\nNEW FILE HAS A DIFFERENT SHAPE")
141 141 print(self.dataShape,newShape,"\n")
142 142 return 0
143 143 else:
144 144 self.dataShape = fp.get('Raw11/Data/Samples/Data').shape[1:]
145 145
146 146
147 147 header = 'Raw11/Data/RadacHeader'
148 148 self.beamCodeByPulse = fp.get(header+'/BeamCode') # LIST OF BEAMS PER PROFILE, TO BE USED ON REARRANGE
149 149 if (self.startDate> datetime.date(2021, 7, 15)): #Se cambió la forma de extracción de Apuntes el 17
150 150 self.beamcodeFile = fp['Setup/Beamcodefile'][()].decode()
151 151 self.trueBeams = self.beamcodeFile.split("\n")
152 152 self.trueBeams.pop()#remove last
153 153 [self.realBeamCode.append(x) for x in self.trueBeams if x not in self.realBeamCode]
154 154 self.beamCode = [int(x, 16) for x in self.realBeamCode]
155 155 else:
156 156 _beamCode= fp.get('Raw11/Data/Beamcodes') #se usa la manera previa al cambio de apuntes
157 157 self.beamCode = _beamCode[0,:]
158 158
159 159
160 160 #self.code = fp.get(header+'/Code') # NOT USE FOR THIS
161 161 self.frameCount = fp.get(header+'/FrameCount')# NOT USE FOR THIS
162 162 self.modeGroup = fp.get(header+'/ModeGroup')# NOT USE FOR THIS
163 163 self.nsamplesPulse = fp.get(header+'/NSamplesPulse')# TO GET NSA OR USING DATA FOR THAT
164 164 self.pulseCount = fp.get(header+'/PulseCount')# NOT USE FOR THIS
165 165 self.radacTime = fp.get(header+'/RadacTime')# 1st TIME ON FILE ANDE CALCULATE THE REST WITH IPP*nindexprofile
166 166 self.timeCount = fp.get(header+'/TimeCount')# NOT USE FOR THIS
167 167 self.timeStatus = fp.get(header+'/TimeStatus')# NOT USE FOR THIS
168 168 self.rangeFromFile = fp.get('Raw11/Data/Samples/Range')
169 169 self.frequency = fp.get('Rx/Frequency')
170 170 txAus = fp.get('Raw11/Data/Pulsewidth')
171 171
172 172
173 173 self.nblocks = self.pulseCount.shape[0] #nblocks
174 174
175 175 self.nprofiles = self.pulseCount.shape[1] #nprofile
176 176 self.nsa = self.nsamplesPulse[0,0] #ngates
177 177 self.nchannels = len(self.beamCode)
178 178 self.ippSeconds = (self.radacTime[0][1] -self.radacTime[0][0]) #Ipp in seconds
179 179 #self.__waitForNewFile = self.nblocks # wait depending on the number of blocks since each block is 1 sec
180 180 self.__waitForNewFile = self.nblocks * self.nprofiles * self.ippSeconds # wait until new file is created
181 181
182 182 #filling radar controller header parameters
183 183 self.__ippKm = self.ippSeconds *.15*1e6 # in km
184 184 self.__txA = (txAus.value)*.15 #(ipp[us]*.15km/1us) in km
185 185 self.__txB = 0
186 186 nWindows=1
187 187 self.__nSamples = self.nsa
188 188 self.__firstHeight = self.rangeFromFile[0][0]/1000 #in km
189 189 self.__deltaHeight = (self.rangeFromFile[0][1] - self.rangeFromFile[0][0])/1000
190 190
191 191 #for now until understand why the code saved is different (code included even though code not in tuf file)
192 192 #self.__codeType = 0
193 193 # self.__nCode = None
194 194 # self.__nBaud = None
195 195 self.__code = self.code
196 196 self.__codeType = 0
197 197 if self.code != None:
198 198 self.__codeType = 1
199 199 self.__nCode = self.nCode
200 200 self.__nBaud = self.nBaud
201 201 #self.__code = 0
202 202
203 203 #filling system header parameters
204 204 self.__nSamples = self.nsa
205 205 self.newProfiles = self.nprofiles/self.nchannels
206 206 self.__channelList = list(range(self.nchannels))
207 207
208 208 self.__frequency = self.frequency[0][0]
209 209
210 210
211 211 return 1
212 212
213 213
214 214 def createBuffers(self):
215 215
216 216 pass
217 217
218 218 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
219 219 self.path = path
220 220 self.startDate = startDate
221 221 self.endDate = endDate
222 222 self.startTime = startTime
223 223 self.endTime = endTime
224 224 self.walk = walk
225 225
226 226 def __checkPath(self):
227 227 if os.path.exists(self.path):
228 228 self.status = 1
229 229 else:
230 230 self.status = 0
231 231 print('Path:%s does not exists'%self.path)
232 232
233 233 return
234 234
235 235
236 236 def __selDates(self, amisr_dirname_format):
237 237 try:
238 238 year = int(amisr_dirname_format[0:4])
239 239 month = int(amisr_dirname_format[4:6])
240 240 dom = int(amisr_dirname_format[6:8])
241 241 thisDate = datetime.date(year,month,dom)
242 242
243 243 if (thisDate>=self.startDate and thisDate <= self.endDate):
244 244 return amisr_dirname_format
245 245 except:
246 246 return None
247 247
248 248
249 249 def __findDataForDates(self,online=False):
250 250
251 251 if not(self.status):
252 252 return None
253 253
254 254 pat = '\d+.\d+'
255 255 dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
256 256 dirnameList = [x for x in dirnameList if x!=None]
257 257 dirnameList = [x.string for x in dirnameList]
258 258 if not(online):
259 259 dirnameList = [self.__selDates(x) for x in dirnameList]
260 260 dirnameList = [x for x in dirnameList if x!=None]
261 261 if len(dirnameList)>0:
262 262 self.status = 1
263 263 self.dirnameList = dirnameList
264 264 self.dirnameList.sort()
265 265 else:
266 266 self.status = 0
267 267 return None
268 268
269 269 def __getTimeFromData(self):
270 270 startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
271 271 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
272 272
273 273 print('Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader))
274 274 print('........................................')
275 275 filter_filenameList = []
276 276 self.filenameList.sort()
277 277 #for i in range(len(self.filenameList)-1):
278 278 for i in range(len(self.filenameList)):
279 279 filename = self.filenameList[i]
280 280 fp = h5py.File(filename,'r')
281 281 time_str = fp.get('Time/RadacTimeString')
282 282
283 283 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
284 284 #startDateTimeStr_File = "2019-12-16 09:21:11"
285 285 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
286 286 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
287 287
288 288 #endDateTimeStr_File = "2019-12-16 11:10:11"
289 289 endDateTimeStr_File = time_str[-1][-1].decode('UTF-8').split('.')[0]
290 290 junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
291 291 endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
292 292
293 293 fp.close()
294 294
295 295 #print("check time", startDateTime_File)
296 296 if self.timezone == 'lt':
297 297 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
298 298 endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
299 299 if (endDateTime_File>=startDateTime_Reader and endDateTime_File<=endDateTime_Reader):
300 300 filter_filenameList.append(filename)
301 301
302 302 if (endDateTime_File>endDateTime_Reader):
303 303 break
304 304
305 305
306 306 filter_filenameList.sort()
307 307 self.filenameList = filter_filenameList
308 308 return 1
309 309
310 310 def __filterByGlob1(self, dirName):
311 311 filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
312 312 filter_files.sort()
313 313 filterDict = {}
314 314 filterDict.setdefault(dirName)
315 315 filterDict[dirName] = filter_files
316 316 return filterDict
317 317
318 318 def __getFilenameList(self, fileListInKeys, dirList):
319 319 for value in fileListInKeys:
320 320 dirName = list(value.keys())[0]
321 321 for file in value[dirName]:
322 322 filename = os.path.join(dirName, file)
323 323 self.filenameList.append(filename)
324 324
325 325
326 326 def __selectDataForTimes(self, online=False):
327 327 #aun no esta implementado el filtro for tiempo
328 328 if not(self.status):
329 329 return None
330 330
331 331 dirList = [os.path.join(self.path,x) for x in self.dirnameList]
332 332
333 333 fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
334 334
335 335 self.__getFilenameList(fileListInKeys, dirList)
336 336 if not(online):
337 337 #filtro por tiempo
338 338 if not(self.all):
339 339 self.__getTimeFromData()
340 340
341 341 if len(self.filenameList)>0:
342 342 self.status = 1
343 343 self.filenameList.sort()
344 344 else:
345 345 self.status = 0
346 346 return None
347 347
348 348 else:
349 349 #get the last file - 1
350 self.filenameList = [self.filenameList[-1]]
350 self.filenameList = [self.filenameList[-2]]
351 351 new_dirnameList = []
352 352 for dirname in self.dirnameList:
353 353 junk = numpy.array([dirname in x for x in self.filenameList])
354 354 junk_sum = junk.sum()
355 355 if junk_sum > 0:
356 356 new_dirnameList.append(dirname)
357 357 self.dirnameList = new_dirnameList
358 358 return 1
359 359
360 360 def searchFilesOnLine(self, path, startDate, endDate, startTime=datetime.time(0,0,0),
361 361 endTime=datetime.time(23,59,59),walk=True):
362 362
363 363 if endDate ==None:
364 364 startDate = datetime.datetime.utcnow().date()
365 365 endDate = datetime.datetime.utcnow().date()
366 366
367 367 self.__setParameters(path=path, startDate=startDate, endDate=endDate,startTime = startTime,endTime=endTime, walk=walk)
368 368
369 369 self.__checkPath()
370 370
371 371 self.__findDataForDates(online=True)
372 372
373 373 self.dirnameList = [self.dirnameList[-1]]
374 374
375 375 self.__selectDataForTimes(online=True)
376 376
377 377 return
378 378
379 379
380 380 def searchFilesOffLine(self,
381 381 path,
382 382 startDate,
383 383 endDate,
384 384 startTime=datetime.time(0,0,0),
385 385 endTime=datetime.time(23,59,59),
386 386 walk=True):
387 387
388 388 self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
389 389
390 390 self.__checkPath()
391 391
392 392 self.__findDataForDates()
393 393
394 394 self.__selectDataForTimes()
395 395
396 396 for i in range(len(self.filenameList)):
397 397 print("%s" %(self.filenameList[i]))
398 398
399 399 return
400 400
401 401 def __setNextFileOffline(self):
402 402
403 403 try:
404 404 self.filename = self.filenameList[self.fileIndex]
405 405 self.amisrFilePointer = h5py.File(self.filename,'r')
406 406 self.fileIndex += 1
407 407 except:
408 408 self.flagNoMoreFiles = 1
409 409 print("No more Files")
410 410 return 0
411 411
412 412 self.flagIsNewFile = 1
413 413 print("Setting the file: %s"%self.filename)
414 414
415 415 return 1
416 416
417 417
418 418 def __setNextFileOnline(self):
419 419 filename = self.filenameList[0]
420 420 if self.__filename_online != None:
421 421 self.__selectDataForTimes(online=True)
422 422 filename = self.filenameList[0]
423 423 wait = 0
424 #self.__waitForNewFile=180 ## DEBUG:
424 self.__waitForNewFile=300 ## DEBUG:
425 425 while self.__filename_online == filename:
426 426 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
427 427 if wait == 5:
428 428 self.flagNoMoreFiles = 1
429 429 return 0
430 430 sleep(self.__waitForNewFile)
431 431 self.__selectDataForTimes(online=True)
432 432 filename = self.filenameList[0]
433 433 wait += 1
434 434
435 435 self.__filename_online = filename
436 436
437 437 self.amisrFilePointer = h5py.File(filename,'r')
438 438 self.flagIsNewFile = 1
439 439 self.filename = filename
440 440 print("Setting the file: %s"%self.filename)
441 441 return 1
442 442
443 443
444 444 def readData(self):
445 445 buffer = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
446 446 re = buffer[:,:,:,0]
447 447 im = buffer[:,:,:,1]
448 448 dataset = re + im*1j
449 449
450 450 self.radacTime = self.amisrFilePointer.get('Raw11/Data/RadacHeader/RadacTime')
451 451 timeset = self.radacTime[:,0]
452 452
453 453 return dataset,timeset
454 454
455 455 def reshapeData(self):
456 456 #self.beamCodeByPulse, self.beamCode, self.nblocks, self.nprofiles, self.nsa,
457 457 channels = self.beamCodeByPulse[0,:]
458 458 nchan = self.nchannels
459 459 #self.newProfiles = self.nprofiles/nchan #must be defined on filljroheader
460 460 nblocks = self.nblocks
461 461 nsamples = self.nsa
462 462
463 463 #Dimensions : nChannels, nProfiles, nSamples
464 464 new_block = numpy.empty((nblocks, nchan, numpy.int_(self.newProfiles), nsamples), dtype="complex64")
465 465 ############################################
466 466
467 467 for thisChannel in range(nchan):
468 468 new_block[:,thisChannel,:,:] = self.dataset[:,numpy.where(channels==self.beamCode[thisChannel])[0],:]
469 469
470 470
471 471 new_block = numpy.transpose(new_block, (1,0,2,3))
472 472 new_block = numpy.reshape(new_block, (nchan,-1, nsamples))
473 473
474 474 return new_block
475 475
476 476 def updateIndexes(self):
477 477
478 478 pass
479 479
480 480 def fillJROHeader(self):
481 481
482 482 #fill radar controller header
483 483 self.dataOut.radarControllerHeaderObj = RadarControllerHeader(ipp=self.__ippKm,
484 484 txA=self.__txA,
485 485 txB=0,
486 486 nWindows=1,
487 487 nHeights=self.__nSamples,
488 488 firstHeight=self.__firstHeight,
489 489 deltaHeight=self.__deltaHeight,
490 490 codeType=self.__codeType,
491 491 nCode=self.__nCode, nBaud=self.__nBaud,
492 492 code = self.__code,
493 493 fClock=1)
494 494
495 495 #fill system header
496 496 self.dataOut.systemHeaderObj = SystemHeader(nSamples=self.__nSamples,
497 497 nProfiles=self.newProfiles,
498 498 nChannels=len(self.__channelList),
499 499 adcResolution=14,
500 500 pciDioBusWidth=32)
501 501
502 502 self.dataOut.type = "Voltage"
503 503
504 504 self.dataOut.data = None
505 505
506 506 self.dataOut.dtype = numpy.dtype([('real','<i8'),('imag','<i8')])
507 507
508 508 # self.dataOut.nChannels = 0
509 509
510 510 # self.dataOut.nHeights = 0
511 511
512 512 self.dataOut.nProfiles = self.newProfiles*self.nblocks
513 513
514 514 #self.dataOut.heightList = self.__firstHeigth + numpy.arange(self.__nSamples, dtype = numpy.float)*self.__deltaHeigth
515 515 ranges = numpy.reshape(self.rangeFromFile.value,(-1))
516 516 self.dataOut.heightList = ranges/1000.0 #km
517 517
518 518
519 519 self.dataOut.channelList = self.__channelList
520 520
521 521 self.dataOut.blocksize = self.dataOut.nChannels * self.dataOut.nHeights
522 522
523 523 # self.dataOut.channelIndexList = None
524 524
525 525 self.dataOut.flagNoData = True
526 526
527 527 #Set to TRUE if the data is discontinuous
528 528 self.dataOut.flagDiscontinuousBlock = False
529 529
530 530 self.dataOut.utctime = None
531 531
532 532 #self.dataOut.timeZone = -5 #self.__timezone/60 #timezone like jroheader, difference in minutes between UTC and localtime
533 533 if self.timezone == 'lt':
534 534 self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
535 535 else:
536 536 self.dataOut.timeZone = 0 #by default time is UTC
537 537
538 538 self.dataOut.dstFlag = 0
539 539
540 540 self.dataOut.errorCount = 0
541 541
542 542 self.dataOut.nCohInt = 1
543 543
544 544 self.dataOut.flagDecodeData = False #asumo que la data esta decodificada
545 545
546 546 self.dataOut.flagDeflipData = False #asumo que la data esta sin flip
547 547
548 548 self.dataOut.flagShiftFFT = False
549 549
550 550 self.dataOut.ippSeconds = self.ippSeconds
551 551
552 552 #Time interval between profiles
553 553 #self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
554 554
555 555 self.dataOut.frequency = self.__frequency
556 556 self.dataOut.realtime = self.online
557 557 pass
558 558
559 559 def readNextFile(self,online=False):
560 560
561 561 if not(online):
562 562 newFile = self.__setNextFileOffline()
563 563 else:
564 564 newFile = self.__setNextFileOnline()
565 565
566 566 if not(newFile):
567 567 self.dataOut.error = True
568 568 return 0
569 569
570 570 if not self.readAMISRHeader(self.amisrFilePointer):
571 571 self.dataOut.error = True
572 572 return 0
573 573
574 574 self.createBuffers()
575 575 self.fillJROHeader()
576 576
577 577 #self.__firstFile = False
578 578
579 579
580 580
581 581 self.dataset,self.timeset = self.readData()
582 582
583 583 if self.endDate!=None:
584 584 endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
585 585 time_str = self.amisrFilePointer.get('Time/RadacTimeString')
586 586 startDateTimeStr_File = time_str[0][0].decode('UTF-8').split('.')[0]
587 587 junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
588 588 startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
589 589 if self.timezone == 'lt':
590 590 startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
591 591 if (startDateTime_File>endDateTime_Reader):
592 592 return 0
593 593
594 594 self.jrodataset = self.reshapeData()
595 595 #----self.updateIndexes()
596 596 self.profileIndex = 0
597 597
598 598 return 1
599 599
600 600
601 601 def __hasNotDataInBuffer(self):
602 602 if self.profileIndex >= (self.newProfiles*self.nblocks):
603 603 return 1
604 604 return 0
605 605
606 606
607 607 def getData(self):
608 608
609 609 if self.flagNoMoreFiles:
610 610 self.dataOut.flagNoData = True
611 611 return 0
612 612
613 613 if self.__hasNotDataInBuffer():
614 614 if not (self.readNextFile(self.online)):
615 615 return 0
616 616
617 617
618 618 if self.dataset is None: # setear esta condicion cuando no hayan datos por leer
619 619 self.dataOut.flagNoData = True
620 620 return 0
621 621
622 622 #self.dataOut.data = numpy.reshape(self.jrodataset[self.profileIndex,:],(1,-1))
623 623
624 624 self.dataOut.data = self.jrodataset[:,self.profileIndex,:]
625 625
626 626 #print("R_t",self.timeset)
627 627
628 628 #self.dataOut.utctime = self.jrotimeset[self.profileIndex]
629 629 #verificar basic header de jro data y ver si es compatible con este valor
630 630 #self.dataOut.utctime = self.timeset + (self.profileIndex * self.ippSeconds * self.nchannels)
631 631 indexprof = numpy.mod(self.profileIndex, self.newProfiles)
632 632 indexblock = self.profileIndex/self.newProfiles
633 633 #print (indexblock, indexprof)
634 634 diffUTC = 1.8e4 #UTC diference from peru in seconds --Joab
635 635 diffUTC = 0
636 636 t_comp = (indexprof * self.ippSeconds * self.nchannels) + diffUTC #
637 637
638 638 #print("utc :",indexblock," __ ",t_comp)
639 639 #print(numpy.shape(self.timeset))
640 640 self.dataOut.utctime = self.timeset[numpy.int_(indexblock)] + t_comp
641 641 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
642 642 #print(self.dataOut.utctime)
643 643 self.dataOut.profileIndex = self.profileIndex
644 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
644 645 self.dataOut.flagNoData = False
645 646 # if indexprof == 0:
646 647 # print self.dataOut.utctime
647 648
648 649 self.profileIndex += 1
649 650
650 return self.dataOut.data
651 #return self.dataOut.data
651 652
652 653
653 654 def run(self, **kwargs):
654 655 '''
655 656 This method will be called many times so here you should put all your code
656 657 '''
657 658 #print("running kamisr")
658 659 if not self.isConfig:
659 660 self.setup(**kwargs)
660 661 self.isConfig = True
661 662
662 663 self.getData()
664 return(self.dataOut.data)
Show file before | Show file after | Hide comments
@@ -1,347 +1,346
1 1 import numpy
2 2
3 3 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
4 4 from schainpy.model.data.jrodata import SpectraHeis
5 5 from schainpy.utils import log
6 6
7 7
8
9 8 class SpectraHeisProc(ProcessingUnit):
10 9
11 10 def __init__(self):#, **kwargs):
12 11
13 12 ProcessingUnit.__init__(self)#, **kwargs)
14 13
15 14 # self.buffer = None
16 15 # self.firstdatatime = None
17 16 # self.profIndex = 0
18 17 self.dataOut = SpectraHeis()
19 18
20 19 def __updateObjFromVoltage(self):
21 20
22 21 self.dataOut.timeZone = self.dataIn.timeZone
23 22 self.dataOut.dstFlag = self.dataIn.dstFlag
24 23 self.dataOut.errorCount = self.dataIn.errorCount
25 24 self.dataOut.useLocalTime = self.dataIn.useLocalTime
26 25
27 26 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()#
28 27 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()#
29 28 self.dataOut.channelList = self.dataIn.channelList
30 29 self.dataOut.heightList = self.dataIn.heightList
31 30 # self.dataOut.dtype = self.dataIn.dtype
32 31 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
33 32 # self.dataOut.nHeights = self.dataIn.nHeights
34 33 # self.dataOut.nChannels = self.dataIn.nChannels
35 34 self.dataOut.nBaud = self.dataIn.nBaud
36 35 self.dataOut.nCode = self.dataIn.nCode
37 36 self.dataOut.code = self.dataIn.code
38 37 # self.dataOut.nProfiles = 1
39 38 self.dataOut.ippFactor = 1
40 39 self.dataOut.noise_estimation = None
41 40 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
42 41 self.dataOut.nFFTPoints = self.dataIn.nHeights
43 42 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
44 43 # self.dataOut.flagNoData = self.dataIn.flagNoData
45 44 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
46 45 self.dataOut.utctime = self.dataIn.utctime
47 46 # self.dataOut.utctime = self.firstdatatime
48 47 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
49 48 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
50 49 # self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
51 50 self.dataOut.nCohInt = self.dataIn.nCohInt
52 51 self.dataOut.nIncohInt = 1
53 52 # self.dataOut.ippSeconds= self.dataIn.ippSeconds
54 53 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
55 54
56 55 # self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nIncohInt
57 56 # self.dataOut.set=self.dataIn.set
58 57 # self.dataOut.deltaHeight=self.dataIn.deltaHeight
59 58
60 59
61 60 def __updateObjFromFits(self):
62 61
63 62 self.dataOut.utctime = self.dataIn.utctime
64 63 # self.dataOut.channelIndexList = self.dataIn.channelIndexList
65 64
66 65 self.dataOut.channelList = self.dataIn.channelList
67 66 self.dataOut.heightList = self.dataIn.heightList
68 67 self.dataOut.data_spc = self.dataIn.data
69 68 self.dataOut.ippSeconds = self.dataIn.ippSeconds
70 69 self.dataOut.nCohInt = self.dataIn.nCohInt
71 70 self.dataOut.nIncohInt = self.dataIn.nIncohInt
72 71 # self.dataOut.timeInterval = self.dataIn.timeInterval
73 72 self.dataOut.timeZone = self.dataIn.timeZone
74 73 self.dataOut.useLocalTime = True
75 74 # self.dataOut.
76 75 # self.dataOut.
77 76
78 77 def __getFft(self):
79 78
80 79 fft_volt = numpy.fft.fft(self.dataIn.data, axis=1)
81 80 fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
82 81 spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))/(self.dataOut.nFFTPoints)
83 82 self.dataOut.data_spc = spc
84 83
85 84 def run(self):
86 85
87 86 self.dataOut.flagNoData = True
88 87
89 88 if self.dataIn.type == "Fits":
90 89 self.__updateObjFromFits()
91 90 self.dataOut.flagNoData = False
92 91 return
93 92
94 93 if self.dataIn.type == "SpectraHeis":
95 94 self.dataOut.copy(self.dataIn)
96 95 return
97 96
98 97 if self.dataIn.type == "Voltage":
99 98 self.__updateObjFromVoltage()
100 99 self.__getFft()
101 100 self.dataOut.flagNoData = False
102 101
103 102 return
104 103
105 104 raise ValueError("The type object %s is not valid"%(self.dataIn.type))
106 105
107 106
108 107 def selectChannels(self, channelList):
109 108
110 109 channelIndexList = []
111 110
112 111 for channel in channelList:
113 112 index = self.dataOut.channelList.index(channel)
114 113 channelIndexList.append(index)
115 114
116 115 self.selectChannelsByIndex(channelIndexList)
117 116
118 117 def selectChannelsByIndex(self, channelIndexList):
119 118 """
120 119 Selecciona un bloque de datos en base a canales segun el channelIndexList
121 120
122 121 Input:
123 122 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
124 123
125 124 Affected:
126 125 self.dataOut.data
127 126 self.dataOut.channelIndexList
128 127 self.dataOut.nChannels
129 128 self.dataOut.m_ProcessingHeader.totalSpectra
130 129 self.dataOut.systemHeaderObj.numChannels
131 130 self.dataOut.m_ProcessingHeader.blockSize
132 131
133 132 Return:
134 133 None
135 134 """
136 135
137 136 for channelIndex in channelIndexList:
138 137 if channelIndex not in self.dataOut.channelIndexList:
139 138 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
140 139
141 140 data_spc = self.dataOut.data_spc[channelIndexList,:]
142 141
143 142 self.dataOut.data_spc = data_spc
144 143 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
145 144
146 145 return 1
147 146
148 147
149 148 class IncohInt4SpectraHeis(Operation):
150 149
151 150 isConfig = False
152 151
153 152 __profIndex = 0
154 153 __withOverapping = False
155 154
156 155 __byTime = False
157 156 __initime = None
158 157 __lastdatatime = None
159 158 __integrationtime = None
160 159
161 160 __buffer = None
162 161
163 162 __dataReady = False
164 163
165 164 n = None
166 165
167 166 def __init__(self):#, **kwargs):
168 167
169 168 Operation.__init__(self)#, **kwargs)
170 169 # self.isConfig = False
171 170
172 171 def setup(self, n=None, timeInterval=None, overlapping=False):
173 172 """
174 173 Set the parameters of the integration class.
175 174
176 175 Inputs:
177 176
178 177 n : Number of coherent integrations
179 178 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
180 179 overlapping :
181 180
182 181 """
183 182
184 183 self.__initime = None
185 184 self.__lastdatatime = 0
186 185 self.__buffer = None
187 186 self.__dataReady = False
188 187
189 188
190 189 if n == None and timeInterval == None:
191 190 raise ValueError("n or timeInterval should be specified ...")
192 191
193 192 if n != None:
194 193 self.n = n
195 194 self.__byTime = False
196 195 else:
197 196 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
198 197 self.n = 9999
199 198 self.__byTime = True
200 199
201 200 if overlapping:
202 201 self.__withOverapping = True
203 202 self.__buffer = None
204 203 else:
205 204 self.__withOverapping = False
206 205 self.__buffer = 0
207 206
208 207 self.__profIndex = 0
209 208
210 209 def putData(self, data):
211 210
212 211 """
213 212 Add a profile to the __buffer and increase in one the __profileIndex
214 213
215 214 """
216 215
217 216 if not self.__withOverapping:
218 217 self.__buffer += data.copy()
219 218 self.__profIndex += 1
220 219 return
221 220
222 221 #Overlapping data
223 222 nChannels, nHeis = data.shape
224 223 data = numpy.reshape(data, (1, nChannels, nHeis))
225 224
226 225 #If the buffer is empty then it takes the data value
227 226 if self.__buffer is None:
228 227 self.__buffer = data
229 228 self.__profIndex += 1
230 229 return
231 230
232 231 #If the buffer length is lower than n then stakcing the data value
233 232 if self.__profIndex < self.n:
234 233 self.__buffer = numpy.vstack((self.__buffer, data))
235 234 self.__profIndex += 1
236 235 return
237 236
238 237 #If the buffer length is equal to n then replacing the last buffer value with the data value
239 238 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
240 239 self.__buffer[self.n-1] = data
241 240 self.__profIndex = self.n
242 241 return
243 242
244 243
245 244 def pushData(self):
246 245 """
247 246 Return the sum of the last profiles and the profiles used in the sum.
248 247
249 248 Affected:
250 249
251 250 self.__profileIndex
252 251
253 252 """
254 253
255 254 if not self.__withOverapping:
256 255 data = self.__buffer
257 256 n = self.__profIndex
258 257
259 258 self.__buffer = 0
260 259 self.__profIndex = 0
261 260
262 261 return data, n
263 262
264 263 #Integration with Overlapping
265 264 data = numpy.sum(self.__buffer, axis=0)
266 265 n = self.__profIndex
267 266
268 267 return data, n
269 268
270 269 def byProfiles(self, data):
271 270
272 271 self.__dataReady = False
273 272 avgdata = None
274 273 # n = None
275 274
276 275 self.putData(data)
277 276
278 277 if self.__profIndex == self.n:
279 278
280 279 avgdata, n = self.pushData()
281 280 self.__dataReady = True
282 281
283 282 return avgdata
284 283
285 284 def byTime(self, data, datatime):
286 285
287 286 self.__dataReady = False
288 287 avgdata = None
289 288 n = None
290 289
291 290 self.putData(data)
292 291
293 292 if (datatime - self.__initime) >= self.__integrationtime:
294 293 avgdata, n = self.pushData()
295 294 self.n = n
296 295 self.__dataReady = True
297 296
298 297 return avgdata
299 298
300 299 def integrate(self, data, datatime=None):
301 300
302 301 if self.__initime == None:
303 302 self.__initime = datatime
304 303
305 304 if self.__byTime:
306 305 avgdata = self.byTime(data, datatime)
307 306 else:
308 307 avgdata = self.byProfiles(data)
309 308
310 309
311 310 self.__lastdatatime = datatime
312 311
313 312 if avgdata is None:
314 313 return None, None
315 314
316 315 avgdatatime = self.__initime
317 316
318 317 deltatime = datatime -self.__lastdatatime
319 318
320 319 if not self.__withOverapping:
321 320 self.__initime = datatime
322 321 else:
323 322 self.__initime += deltatime
324 323
325 324 return avgdata, avgdatatime
326 325
327 326 def run(self, dataOut, n=None, timeInterval=None, overlapping=False, **kwargs):
328 327
329 328 if not self.isConfig:
330 329 self.setup(n=n, timeInterval=timeInterval, overlapping=overlapping)
331 330 self.isConfig = True
332 331
333 332 avgdata, avgdatatime = self.integrate(dataOut.data_spc, dataOut.utctime)
334 333
335 334 # dataOut.timeInterval *= n
336 335 dataOut.flagNoData = True
337 336
338 337 if self.__dataReady:
339 338 dataOut.data_spc = avgdata
340 339 dataOut.nIncohInt *= self.n
341 340 # dataOut.nCohInt *= self.n
342 341 dataOut.utctime = avgdatatime
343 342 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
344 343 # dataOut.timeInterval = self.__timeInterval*self.n
345 344 dataOut.flagNoData = False
346 345
347 346 return dataOut
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (559 lines changed) Show them Hide them
@@ -1,898 +1,1457
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 import math
15 16
16 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 18 from schainpy.model.data.jrodata import Spectra
18 19 from schainpy.model.data.jrodata import hildebrand_sekhon
19 20 from schainpy.utils import log
20 21
22 from scipy.optimize import curve_fit
23
21 24
22 25 class SpectraProc(ProcessingUnit):
23 26
24 27 def __init__(self):
25 28
26 29 ProcessingUnit.__init__(self)
27 30
28 31 self.buffer = None
29 32 self.firstdatatime = None
30 33 self.profIndex = 0
31 34 self.dataOut = Spectra()
32 35 self.id_min = None
33 36 self.id_max = None
34 37 self.setupReq = False #Agregar a todas las unidades de proc
35 38
36 39 def __updateSpecFromVoltage(self):
37 40
38 41 self.dataOut.timeZone = self.dataIn.timeZone
39 42 self.dataOut.dstFlag = self.dataIn.dstFlag
40 43 self.dataOut.errorCount = self.dataIn.errorCount
41 44 self.dataOut.useLocalTime = self.dataIn.useLocalTime
42 45 try:
43 46 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
44 47 except:
45 48 pass
46 49 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
47 50 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
48 51 self.dataOut.channelList = self.dataIn.channelList
49 52 self.dataOut.heightList = self.dataIn.heightList
50 53 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
51 54 self.dataOut.nProfiles = self.dataOut.nFFTPoints
52 55 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
53 56 self.dataOut.utctime = self.firstdatatime
54 57 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
55 58 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
56 59 self.dataOut.flagShiftFFT = False
57 60 self.dataOut.nCohInt = self.dataIn.nCohInt
58 61 self.dataOut.nIncohInt = 1
59 62 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
60 63 self.dataOut.frequency = self.dataIn.frequency
61 64 self.dataOut.realtime = self.dataIn.realtime
62 65 self.dataOut.azimuth = self.dataIn.azimuth
63 66 self.dataOut.zenith = self.dataIn.zenith
64 67 self.dataOut.beam.codeList = self.dataIn.beam.codeList
65 68 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
66 69 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
67 70
68 71 def __getFft(self):
69 72 """
70 73 Convierte valores de Voltaje a Spectra
71 74
72 75 Affected:
73 76 self.dataOut.data_spc
74 77 self.dataOut.data_cspc
75 78 self.dataOut.data_dc
76 79 self.dataOut.heightList
77 80 self.profIndex
78 81 self.buffer
79 82 self.dataOut.flagNoData
80 83 """
81 84 fft_volt = numpy.fft.fft(
82 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
83 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
84 87 dc = fft_volt[:, 0, :]
85 88
86 89 # calculo de self-spectra
87 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
88 91 spc = fft_volt * numpy.conjugate(fft_volt)
89 92 spc = spc.real
90 93
91 94 blocksize = 0
92 95 blocksize += dc.size
93 96 blocksize += spc.size
94 97
95 98 cspc = None
96 99 pairIndex = 0
97 100 if self.dataOut.pairsList != None:
98 101 # calculo de cross-spectra
99 102 cspc = numpy.zeros(
100 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
101 104 for pair in self.dataOut.pairsList:
102 105 if pair[0] not in self.dataOut.channelList:
103 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
104 107 str(pair), str(self.dataOut.channelList)))
105 108 if pair[1] not in self.dataOut.channelList:
106 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
107 110 str(pair), str(self.dataOut.channelList)))
108 111
109 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
110 113 numpy.conjugate(fft_volt[pair[1], :, :])
111 114 pairIndex += 1
112 115 blocksize += cspc.size
113 116
114 117 self.dataOut.data_spc = spc
115 118 self.dataOut.data_cspc = cspc
116 119 self.dataOut.data_dc = dc
117 120 self.dataOut.blockSize = blocksize
118 121 self.dataOut.flagShiftFFT = False
119 122
120 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
121 124
122 125 if self.dataIn.type == "Spectra":
123 126 self.dataOut.copy(self.dataIn)
124 127 if shift_fft:
125 128 #desplaza a la derecha en el eje 2 determinadas posiciones
126 129 shift = int(self.dataOut.nFFTPoints/2)
127 130 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
128 131
129 132 if self.dataOut.data_cspc is not None:
130 133 #desplaza a la derecha en el eje 2 determinadas posiciones
131 134 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
132 135 if pairsList:
133 136 self.__selectPairs(pairsList)
134 137
135 138 elif self.dataIn.type == "Voltage":
136 139
137 140 self.dataOut.flagNoData = True
138 141
139 142 if nFFTPoints == None:
140 143 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
141 144
142 145 if nProfiles == None:
143 146 nProfiles = nFFTPoints
144 147
145 148 if ippFactor == None:
146 149 self.dataOut.ippFactor = 1
147 150
148 151 self.dataOut.nFFTPoints = nFFTPoints
149 152
150 153 if self.buffer is None:
151 154 self.buffer = numpy.zeros((self.dataIn.nChannels,
152 155 nProfiles,
153 156 self.dataIn.nHeights),
154 157 dtype='complex')
155 158
156 159 if self.dataIn.flagDataAsBlock:
157 160 nVoltProfiles = self.dataIn.data.shape[1]
158 161
159 162 if nVoltProfiles == nProfiles:
160 163 self.buffer = self.dataIn.data.copy()
161 164 self.profIndex = nVoltProfiles
162 165
163 166 elif nVoltProfiles < nProfiles:
164 167
165 168 if self.profIndex == 0:
166 169 self.id_min = 0
167 170 self.id_max = nVoltProfiles
168 171
169 172 self.buffer[:, self.id_min:self.id_max,
170 173 :] = self.dataIn.data
171 174 self.profIndex += nVoltProfiles
172 175 self.id_min += nVoltProfiles
173 176 self.id_max += nVoltProfiles
174 177 else:
175 178 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
176 179 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
177 180 self.dataOut.flagNoData = True
178 181 else:
179 182 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
180 183 self.profIndex += 1
181 184
182 185 if self.firstdatatime == None:
183 186 self.firstdatatime = self.dataIn.utctime
184 187
185 188 if self.profIndex == nProfiles:
186 189 self.__updateSpecFromVoltage()
187 190 if pairsList == None:
188 191 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
189 192 else:
190 193 self.dataOut.pairsList = pairsList
191 194 self.__getFft()
192 195 self.dataOut.flagNoData = False
193 196 self.firstdatatime = None
194 197 self.profIndex = 0
195 198 else:
196 199 raise ValueError("The type of input object '%s' is not valid".format(
197 200 self.dataIn.type))
198 201
199 202 def __selectPairs(self, pairsList):
200 203
201 204 if not pairsList:
202 205 return
203 206
204 207 pairs = []
205 208 pairsIndex = []
206 209
207 210 for pair in pairsList:
208 211 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
209 212 continue
210 213 pairs.append(pair)
211 214 pairsIndex.append(pairs.index(pair))
212 215
213 216 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
214 217 self.dataOut.pairsList = pairs
215 218
216 219 return
217 220
218 221 def selectFFTs(self, minFFT, maxFFT ):
219 222 """
220 223 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
221 224 minFFT<= FFT <= maxFFT
222 225 """
223 226
224 227 if (minFFT > maxFFT):
225 228 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
226 229
227 230 if (minFFT < self.dataOut.getFreqRange()[0]):
228 231 minFFT = self.dataOut.getFreqRange()[0]
229 232
230 233 if (maxFFT > self.dataOut.getFreqRange()[-1]):
231 234 maxFFT = self.dataOut.getFreqRange()[-1]
232 235
233 236 minIndex = 0
234 237 maxIndex = 0
235 238 FFTs = self.dataOut.getFreqRange()
236 239
237 240 inda = numpy.where(FFTs >= minFFT)
238 241 indb = numpy.where(FFTs <= maxFFT)
239 242
240 243 try:
241 244 minIndex = inda[0][0]
242 245 except:
243 246 minIndex = 0
244 247
245 248 try:
246 249 maxIndex = indb[0][-1]
247 250 except:
248 251 maxIndex = len(FFTs)
249 252
250 253 self.selectFFTsByIndex(minIndex, maxIndex)
251 254
252 255 return 1
253 256
254 257 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
255 258 newheis = numpy.where(
256 259 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
257 260
258 261 if hei_ref != None:
259 262 newheis = numpy.where(self.dataOut.heightList > hei_ref)
260 263
261 264 minIndex = min(newheis[0])
262 265 maxIndex = max(newheis[0])
263 266 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
264 267 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
265 268
266 269 # determina indices
267 270 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
268 271 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
269 272 avg_dB = 10 * \
270 273 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
271 274 beacon_dB = numpy.sort(avg_dB)[-nheis:]
272 275 beacon_heiIndexList = []
273 276 for val in avg_dB.tolist():
274 277 if val >= beacon_dB[0]:
275 278 beacon_heiIndexList.append(avg_dB.tolist().index(val))
276 279
277 280 #data_spc = data_spc[:,:,beacon_heiIndexList]
278 281 data_cspc = None
279 282 if self.dataOut.data_cspc is not None:
280 283 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
281 284 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
282 285
283 286 data_dc = None
284 287 if self.dataOut.data_dc is not None:
285 288 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
286 289 #data_dc = data_dc[:,beacon_heiIndexList]
287 290
288 291 self.dataOut.data_spc = data_spc
289 292 self.dataOut.data_cspc = data_cspc
290 293 self.dataOut.data_dc = data_dc
291 294 self.dataOut.heightList = heightList
292 295 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
293 296
294 297 return 1
295 298
296 299 def selectFFTsByIndex(self, minIndex, maxIndex):
297 300 """
298 301
299 302 """
300 303
301 304 if (minIndex < 0) or (minIndex > maxIndex):
302 305 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
303 306
304 307 if (maxIndex >= self.dataOut.nProfiles):
305 308 maxIndex = self.dataOut.nProfiles-1
306 309
307 310 #Spectra
308 311 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
309 312
310 313 data_cspc = None
311 314 if self.dataOut.data_cspc is not None:
312 315 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
313 316
314 317 data_dc = None
315 318 if self.dataOut.data_dc is not None:
316 319 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
317 320
318 321 self.dataOut.data_spc = data_spc
319 322 self.dataOut.data_cspc = data_cspc
320 323 self.dataOut.data_dc = data_dc
321 324
322 325 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
323 326 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
324 327 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
325 328
326 329 return 1
327 330
328 331 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
329 332 # validacion de rango
330 333 if minHei == None:
331 334 minHei = self.dataOut.heightList[0]
332 335
333 336 if maxHei == None:
334 337 maxHei = self.dataOut.heightList[-1]
335 338
336 339 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
337 340 print('minHei: %.2f is out of the heights range' % (minHei))
338 341 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
339 342 minHei = self.dataOut.heightList[0]
340 343
341 344 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
342 345 print('maxHei: %.2f is out of the heights range' % (maxHei))
343 346 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
344 347 maxHei = self.dataOut.heightList[-1]
345 348
346 349 # validacion de velocidades
347 350 velrange = self.dataOut.getVelRange(1)
348 351
349 352 if minVel == None:
350 353 minVel = velrange[0]
351 354
352 355 if maxVel == None:
353 356 maxVel = velrange[-1]
354 357
355 358 if (minVel < velrange[0]) or (minVel > maxVel):
356 359 print('minVel: %.2f is out of the velocity range' % (minVel))
357 360 print('minVel is setting to %.2f' % (velrange[0]))
358 361 minVel = velrange[0]
359 362
360 363 if (maxVel > velrange[-1]) or (maxVel < minVel):
361 364 print('maxVel: %.2f is out of the velocity range' % (maxVel))
362 365 print('maxVel is setting to %.2f' % (velrange[-1]))
363 366 maxVel = velrange[-1]
364 367
365 368 # seleccion de indices para rango
366 369 minIndex = 0
367 370 maxIndex = 0
368 371 heights = self.dataOut.heightList
369 372
370 373 inda = numpy.where(heights >= minHei)
371 374 indb = numpy.where(heights <= maxHei)
372 375
373 376 try:
374 377 minIndex = inda[0][0]
375 378 except:
376 379 minIndex = 0
377 380
378 381 try:
379 382 maxIndex = indb[0][-1]
380 383 except:
381 384 maxIndex = len(heights)
382 385
383 386 if (minIndex < 0) or (minIndex > maxIndex):
384 387 raise ValueError("some value in (%d,%d) is not valid" % (
385 388 minIndex, maxIndex))
386 389
387 390 if (maxIndex >= self.dataOut.nHeights):
388 391 maxIndex = self.dataOut.nHeights - 1
389 392
390 393 # seleccion de indices para velocidades
391 394 indminvel = numpy.where(velrange >= minVel)
392 395 indmaxvel = numpy.where(velrange <= maxVel)
393 396 try:
394 397 minIndexVel = indminvel[0][0]
395 398 except:
396 399 minIndexVel = 0
397 400
398 401 try:
399 402 maxIndexVel = indmaxvel[0][-1]
400 403 except:
401 404 maxIndexVel = len(velrange)
402 405
403 406 # seleccion del espectro
404 407 data_spc = self.dataOut.data_spc[:,
405 408 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
406 409 # estimacion de ruido
407 410 noise = numpy.zeros(self.dataOut.nChannels)
408 411
409 412 for channel in range(self.dataOut.nChannels):
410 413 daux = data_spc[channel, :, :]
411 414 sortdata = numpy.sort(daux, axis=None)
412 415 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
413 416
414 417 self.dataOut.noise_estimation = noise.copy()
415 418
416 419 return 1
417 420
418 421 class removeDC(Operation):
419 422
420 423 def run(self, dataOut, mode=2):
421 424 self.dataOut = dataOut
422 425 jspectra = self.dataOut.data_spc
423 426 jcspectra = self.dataOut.data_cspc
424 427
425 428 num_chan = jspectra.shape[0]
426 429 num_hei = jspectra.shape[2]
427 430
428 431 if jcspectra is not None:
429 432 jcspectraExist = True
430 433 num_pairs = jcspectra.shape[0]
431 434 else:
432 435 jcspectraExist = False
433 436
434 437 freq_dc = int(jspectra.shape[1] / 2)
435 438 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
436 439 ind_vel = ind_vel.astype(int)
437 440
438 441 if ind_vel[0] < 0:
439 442 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
440 443
441 444 if mode == 1:
442 445 jspectra[:, freq_dc, :] = (
443 446 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
444 447
445 448 if jcspectraExist:
446 449 jcspectra[:, freq_dc, :] = (
447 450 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
448 451
449 452 if mode == 2:
450 453
451 454 vel = numpy.array([-2, -1, 1, 2])
452 455 xx = numpy.zeros([4, 4])
453 456
454 457 for fil in range(4):
455 458 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
456 459
457 460 xx_inv = numpy.linalg.inv(xx)
458 461 xx_aux = xx_inv[0, :]
459 462
460 463 for ich in range(num_chan):
461 464 yy = jspectra[ich, ind_vel, :]
462 465 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
463 466
464 467 junkid = jspectra[ich, freq_dc, :] <= 0
465 468 cjunkid = sum(junkid)
466 469
467 470 if cjunkid.any():
468 471 jspectra[ich, freq_dc, junkid.nonzero()] = (
469 472 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
470 473
471 474 if jcspectraExist:
472 475 for ip in range(num_pairs):
473 476 yy = jcspectra[ip, ind_vel, :]
474 477 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
475 478
476 479 self.dataOut.data_spc = jspectra
477 480 self.dataOut.data_cspc = jcspectra
478 481
479 482 return self.dataOut
480 483
484 # import matplotlib.pyplot as plt
485
486 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
487 z = (x - a1) / a2
488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
489 return y
490 class CleanRayleigh(Operation):
491
492 def __init__(self):
493
494 Operation.__init__(self)
495 self.i=0
496 self.isConfig = False
497 self.__dataReady = False
498 self.__profIndex = 0
499 self.byTime = False
500 self.byProfiles = False
501
502 self.bloques = None
503 self.bloque0 = None
504
505 self.index = 0
506
507 self.buffer = 0
508 self.buffer2 = 0
509 self.buffer3 = 0
510 #self.min_hei = None
511 #self.max_hei = None
512
513 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
514
515 self.nChannels = dataOut.nChannels
516 self.nProf = dataOut.nProfiles
517 self.nPairs = dataOut.data_cspc.shape[0]
518 self.pairsArray = numpy.array(dataOut.pairsList)
519 self.spectra = dataOut.data_spc
520 self.cspectra = dataOut.data_cspc
521 self.heights = dataOut.heightList #alturas totales
522 self.nHeights = len(self.heights)
523 self.min_hei = min_hei
524 self.max_hei = max_hei
525 if (self.min_hei == None):
526 self.min_hei = 0
527 if (self.max_hei == None):
528 self.max_hei = dataOut.heightList[-1]
529 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
530 self.heightsClean = self.heights[self.hval] #alturas filtradas
531 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
532 self.nHeightsClean = len(self.heightsClean)
533 self.channels = dataOut.channelList
534 self.nChan = len(self.channels)
535 self.nIncohInt = dataOut.nIncohInt
536 self.__initime = dataOut.utctime
537 self.maxAltInd = self.hval[-1]+1
538 self.minAltInd = self.hval[0]
539
540 self.crosspairs = dataOut.pairsList
541 self.nPairs = len(self.crosspairs)
542 self.normFactor = dataOut.normFactor
543 self.nFFTPoints = dataOut.nFFTPoints
544 self.ippSeconds = dataOut.ippSeconds
545 self.currentTime = self.__initime
546 self.pairsArray = numpy.array(dataOut.pairsList)
547 self.factor_stdv = factor_stdv
548
549
550
551 if n != None :
552 self.byProfiles = True
553 self.nIntProfiles = n
554 else:
555 self.__integrationtime = timeInterval
556
557 self.__dataReady = False
558 self.isConfig = True
559
560
561
562 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
563 #print (dataOut.utctime)
564 if not self.isConfig :
565 #print("Setting config")
566 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
567 #print("Config Done")
568 tini=dataOut.utctime
569
570 if self.byProfiles:
571 if self.__profIndex == self.nIntProfiles:
572 self.__dataReady = True
573 else:
574 if (tini - self.__initime) >= self.__integrationtime:
575 #print(tini - self.__initime,self.__profIndex)
576 self.__dataReady = True
577 self.__initime = tini
578
579 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
580
581 if self.__dataReady:
582 print("Data ready",self.__profIndex)
583 self.__profIndex = 0
584 jspc = self.buffer
585 jcspc = self.buffer2
586 #jnoise = self.buffer3
587 self.buffer = dataOut.data_spc
588 self.buffer2 = dataOut.data_cspc
589 #self.buffer3 = dataOut.noise
590 self.currentTime = dataOut.utctime
591 if numpy.any(jspc) :
592 #print( jspc.shape, jcspc.shape)
593 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
594 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
595 self.__dataReady = False
596 #print( jspc.shape, jcspc.shape)
597 dataOut.flagNoData = False
598 else:
599 dataOut.flagNoData = True
600 self.__dataReady = False
601 return dataOut
602 else:
603 #print( len(self.buffer))
604 if numpy.any(self.buffer):
605 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
606 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
607 self.buffer3 += dataOut.data_dc
608 else:
609 self.buffer = dataOut.data_spc
610 self.buffer2 = dataOut.data_cspc
611 self.buffer3 = dataOut.data_dc
612 #print self.index, self.fint
613 #print self.buffer2.shape
614 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
615 self.__profIndex += 1
616 return dataOut ## NOTE: REV
617
618 # if self.index == 0 and self.fint == 1 :
619 # if jspc != None:
620 # print len(jspc), jspc.shape
621 # jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4))
622 # print jspc.shape
623 # dataOut.flagNoData = True
624 # return dataOut
625 # if path != None:
626 # sys.path.append(path)
627 # self.library = importlib.import_module(file)
628 #
629 # To be inserted as a parameter
630
631 #Set constants
632 #constants = self.library.setConstants(dataOut)
633 #dataOut.constants = constants
634
635 #snrth= 20
636
637 #crosspairs = dataOut.groupList
638 #noise = dataOut.noise
639 #print( nProf,heights)
640 #print( jspc.shape, jspc.shape[0])
641 #print noise
642 #print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise)
643 #jnoise = jnoise/N
644 #noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
645 #print( noise)
646 #power = numpy.sum(spectra, axis=1)
647 #print power[0,:]
648 #print("CROSSPAIRS",crosspairs)
649 #nPairs = len(crosspairs)
650 #print(numpy.shape(dataOut.data_spc))
651
652 #absc = dataOut.abscissaList[:-1]
653
654 #print absc.shape
655 #nBlocks=149
656 #print('spectra', spectra.shape)
657 #print('noise print', crosspairs)
658 #print('spectra', spectra.shape)
659 #print('cspectra', cspectra.shape)
660 #print numpy.array(dataOut.data_pre[1]).shape
661 #spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs)
662
663
664 #index = tini.tm_hour*12+tini.tm_min/5
665 # jspc = jspc/self.nFFTPoints/self.normFactor
666 # jcspc = jcspc/self.nFFTPoints/self.normFactor
667
668
669
670 #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
671 #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei)
672 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
673 #jspectra = tmp_spectra*len(jspc[:,0,0,0])
674 #jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
675
676 dataOut.data_spc = tmp_spectra
677 dataOut.data_cspc = tmp_cspectra
678 dataOut.data_dc = self.buffer3
679 dataOut.nIncohInt *= self.nIntProfiles
680 dataOut.utctime = self.currentTime #tiempo promediado
681 print("Time: ",time.localtime(dataOut.utctime))
682 # dataOut.data_spc = sat_spectra
683 # dataOut.data_cspc = sat_cspectra
684 self.buffer = 0
685 self.buffer2 = 0
686 self.buffer3 = 0
687
688 return dataOut
689
690 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
691 print("OP cleanRayleigh")
692 #import matplotlib.pyplot as plt
693 #for k in range(149):
694
695 rfunc = cspectra.copy() #self.bloques
696 val_spc = spectra*0.0 #self.bloque0*0.0
697 val_cspc = cspectra*0.0 #self.bloques*0.0
698 in_sat_spectra = spectra.copy() #self.bloque0
699 in_sat_cspectra = cspectra.copy() #self.bloques
700
701 raxs = math.ceil(math.sqrt(self.nPairs))
702 caxs = math.ceil(self.nPairs/raxs)
703
704 #print(self.hval)
705 #print numpy.absolute(rfunc[:,0,0,14])
706 for ih in range(self.minAltInd,self.maxAltInd):
707 for ifreq in range(self.nFFTPoints):
708 # fig, axs = plt.subplots(raxs, caxs)
709 # fig2, axs2 = plt.subplots(raxs, caxs)
710 col_ax = 0
711 row_ax = 0
712 for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
713 #print("ii: ",ii)
714 if (col_ax%caxs==0 and col_ax!=0):
715 col_ax = 0
716 row_ax += 1
717 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
718 #print(func2clean.shape)
719 val = (numpy.isfinite(func2clean)==True).nonzero()
720
721 if len(val)>0: #limitador
722 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
723 if min_val <= -40 :
724 min_val = -40
725 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
726 if max_val >= 200 :
727 max_val = 200
728 #print min_val, max_val
729 step = 1
730 #print("Getting bins and the histogram")
731 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
732 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
733 #print(len(y_dist),len(binstep[:-1]))
734 #print(row_ax,col_ax, " ..")
735 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
736 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
737 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
738 parg = [numpy.amax(y_dist),mean,sigma]
739 gauss_fit, covariance = None, None
740 newY = None
741 try :
742 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
743 mode = gauss_fit[1]
744 stdv = gauss_fit[2]
745 #print(" FIT OK",gauss_fit)
746 '''
747 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
748 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
749 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
750 axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
751 except:
752 mode = mean
753 stdv = sigma
754 #print("FIT FAIL")
755
756
757 #print(mode,stdv)
758 #Removing echoes greater than mode + 3*stdv
759 #factor_stdv = 2
760 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
761 #noval tiene los indices que se van a remover
762 #print("Pair ",ii," novals: ",len(noval[0]))
763 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
764 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
765 #print(novall)
766 #print(" ")
767 cross_pairs = self.pairsArray[ii]
768 #Getting coherent echoes which are removed.
769 if len(novall[0]) > 0:
770
771 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
772 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
773 val_cspc[novall[0],ii,ifreq,ih] = 1
774 #print("OUT NOVALL 1")
775 #Removing coherent from ISR data
776
777 #print(spectra[:,ii,ifreq,ih])
778 new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
779 mean_cspc = numpy.mean(new_a)
780 new_b = numpy.delete(spectra[:,cross_pairs[0],ifreq,ih], noval[0])
781 mean_spc0 = numpy.mean(new_b)
782 new_c = numpy.delete(spectra[:,cross_pairs[1],ifreq,ih], noval[0])
783 mean_spc1 = numpy.mean(new_c)
784 spectra[noval,cross_pairs[0],ifreq,ih] = mean_spc0
785 spectra[noval,cross_pairs[1],ifreq,ih] = mean_spc1
786 cspectra[noval,ii,ifreq,ih] = mean_cspc
787
788 '''
789 func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
790 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
791 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
792 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
793 axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
794 '''
795
796 col_ax += 1 #contador de ploteo columnas
797 ##print(col_ax)
798 '''
799 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
800 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
801 fig.suptitle(title)
802 fig2.suptitle(title2)
803 plt.show()'''
804
805 ''' channels = channels
806 cross_pairs = cross_pairs
807 #print("OUT NOVALL 2")
808
809 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
810 vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
811 vcross = numpy.concatenate((vcross0,vcross1),axis=None)
812 #print('vcros =', vcross)
813
814 #Getting coherent echoes which are removed.
815 if len(novall) > 0:
816 #val_spc[novall,ii,ifreq,ih] = 1
817 val_spc[ii,ifreq,ih,novall] = 1
818 if len(vcross) > 0:
819 val_cspc[vcross,ifreq,ih,novall] = 1
820
821 #Removing coherent from ISR data.
822 self.bloque0[ii,ifreq,ih,noval] = numpy.nan
823 if len(vcross) > 0:
824 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
825 '''
826
827 print("Getting average of the spectra and cross-spectra from incoherent echoes.")
828 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
829 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
830 for ih in range(self.nHeights):
831 for ifreq in range(self.nFFTPoints):
832 for ich in range(self.nChan):
833 tmp = spectra[:,ich,ifreq,ih]
834 valid = (numpy.isfinite(tmp[:])==True).nonzero()
835 # if ich == 0 and ifreq == 0 and ih == 17 :
836 # print tmp
837 # print valid
838 # print len(valid[0])
839 #print('TMP',tmp)
840 if len(valid[0]) >0 :
841 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
842 #for icr in range(nPairs):
843 for icr in range(self.nPairs):
844 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
845 valid = (numpy.isfinite(tmp)==True).nonzero()
846 if len(valid[0]) > 0:
847 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
848 '''
849 # print('##########################################################')
850 print("Removing fake coherent echoes (at least 4 points around the point)")
851
852 val_spectra = numpy.sum(val_spc,0)
853 val_cspectra = numpy.sum(val_cspc,0)
854
855 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
856 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
857
858 for i in range(nChan):
859 for j in range(nProf):
860 for k in range(nHeights):
861 if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
862 val_spc[:,i,j,k] = 0.0
863 for i in range(nPairs):
864 for j in range(nProf):
865 for k in range(nHeights):
866 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
867 val_cspc[:,i,j,k] = 0.0
868
869 # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan))
870 # if numpy.isfinite(val_spectra)==str(True):
871 # noval = (val_spectra<1).nonzero()
872 # if len(noval) > 0:
873 # val_spc[:,noval] = 0.0
874 # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights))
875
876 #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf))
877 #if numpy.isfinite(val_cspectra)==str(True):
878 # noval = (val_cspectra<1).nonzero()
879 # if len(noval) > 0:
880 # val_cspc[:,noval] = 0.0
881 # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights))
882 tmp_sat_spectra = spectra.copy()
883 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
884 tmp_sat_cspectra = cspectra.copy()
885 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
886 '''
887 # fig = plt.figure(figsize=(6,5))
888 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
889 # ax = fig.add_axes([left, bottom, width, height])
890 # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
891 # ax.clabel(cp, inline=True,fontsize=10)
892 # plt.show()
893 '''
894 val = (val_spc > 0).nonzero()
895 if len(val[0]) > 0:
896 tmp_sat_spectra[val] = in_sat_spectra[val]
897 val = (val_cspc > 0).nonzero()
898 if len(val[0]) > 0:
899 tmp_sat_cspectra[val] = in_sat_cspectra[val]
900
901 print("Getting average of the spectra and cross-spectra from incoherent echoes 2")
902 sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float)
903 sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex)
904 for ih in range(nHeights):
905 for ifreq in range(nProf):
906 for ich in range(nChan):
907 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
908 valid = (numpy.isfinite(tmp)).nonzero()
909 if len(valid[0]) > 0:
910 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
911
912 for icr in range(nPairs):
913 tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
914 valid = (numpy.isfinite(tmp)).nonzero()
915 if len(valid[0]) > 0:
916 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
917 '''
918 #self.__dataReady= True
919 #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
920 #if not self.__dataReady:
921 #return None, None
922 #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra
923 return out_spectra, out_cspectra
924
925 def REM_ISOLATED_POINTS(self,array,rth):
926 # import matplotlib.pyplot as plt
927 if rth == None :
928 rth = 4
929 print("REM ISO")
930 num_prof = len(array[0,:,0])
931 num_hei = len(array[0,0,:])
932 n2d = len(array[:,0,0])
933
934 for ii in range(n2d) :
935 #print ii,n2d
936 tmp = array[ii,:,:]
937 #print tmp.shape, array[ii,101,:],array[ii,102,:]
938
939 # fig = plt.figure(figsize=(6,5))
940 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
941 # ax = fig.add_axes([left, bottom, width, height])
942 # x = range(num_prof)
943 # y = range(num_hei)
944 # cp = ax.contour(y,x,tmp)
945 # ax.clabel(cp, inline=True,fontsize=10)
946 # plt.show()
947
948 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
949 tmp = numpy.reshape(tmp,num_prof*num_hei)
950 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
951 indxs2 = (tmp > 0).nonzero()
952
953 indxs1 = (indxs1[0])
954 indxs2 = indxs2[0]
955 #indxs1 = numpy.array(indxs1[0])
956 #indxs2 = numpy.array(indxs2[0])
957 indxs = None
958 #print indxs1 , indxs2
959 for iv in range(len(indxs2)):
960 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
961 #print len(indxs2), indv
962 if len(indv[0]) > 0 :
963 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
964 # print indxs
965 indxs = indxs[1:]
966 #print(indxs, len(indxs))
967 if len(indxs) < 4 :
968 array[ii,:,:] = 0.
969 return
970
971 xpos = numpy.mod(indxs ,num_hei)
972 ypos = (indxs / num_hei)
973 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
974 #print sx
975 xpos = xpos[sx]
976 ypos = ypos[sx]
977
978 # *********************************** Cleaning isolated points **********************************
979 ic = 0
980 while True :
981 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
982 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
983 #plt.plot(r)
984 #plt.show()
985 no_coh1 = (numpy.isfinite(r)==True).nonzero()
986 no_coh2 = (r <= rth).nonzero()
987 #print r, no_coh1, no_coh2
988 no_coh1 = numpy.array(no_coh1[0])
989 no_coh2 = numpy.array(no_coh2[0])
990 no_coh = None
991 #print valid1 , valid2
992 for iv in range(len(no_coh2)):
993 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
994 if len(indv[0]) > 0 :
995 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
996 no_coh = no_coh[1:]
997 #print len(no_coh), no_coh
998 if len(no_coh) < 4 :
999 #print xpos[ic], ypos[ic], ic
1000 # plt.plot(r)
1001 # plt.show()
1002 xpos[ic] = numpy.nan
1003 ypos[ic] = numpy.nan
1004
1005 ic = ic + 1
1006 if (ic == len(indxs)) :
1007 break
1008 #print( xpos, ypos)
1009
1010 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1011 #print indxs[0]
1012 if len(indxs[0]) < 4 :
1013 array[ii,:,:] = 0.
1014 return
1015
1016 xpos = xpos[indxs[0]]
1017 ypos = ypos[indxs[0]]
1018 for i in range(0,len(ypos)):
1019 ypos[i]=int(ypos[i])
1020 junk = tmp
1021 tmp = junk*0.0
1022
1023 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1024 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1025
1026 #print array.shape
1027 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1028 #print tmp.shape
1029
1030 # fig = plt.figure(figsize=(6,5))
1031 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1032 # ax = fig.add_axes([left, bottom, width, height])
1033 # x = range(num_prof)
1034 # y = range(num_hei)
1035 # cp = ax.contour(y,x,array[ii,:,:])
1036 # ax.clabel(cp, inline=True,fontsize=10)
1037 # plt.show()
1038 return array
1039
481 1040 class removeInterference(Operation):
482 1041
483 1042 def removeInterference2(self):
484 1043
485 1044 cspc = self.dataOut.data_cspc
486 1045 spc = self.dataOut.data_spc
487 1046 Heights = numpy.arange(cspc.shape[2])
488 1047 realCspc = numpy.abs(cspc)
489 1048
490 1049 for i in range(cspc.shape[0]):
491 1050 LinePower= numpy.sum(realCspc[i], axis=0)
492 1051 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
493 1052 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
494 1053 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
495 1054 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
496 1055 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
497 1056
498 1057
499 1058 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
500 1059 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
501 1060 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
502 1061 cspc[i,InterferenceRange,:] = numpy.NaN
503 1062
504 1063 self.dataOut.data_cspc = cspc
505 1064
506 1065 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
507 1066
508 1067 jspectra = self.dataOut.data_spc
509 1068 jcspectra = self.dataOut.data_cspc
510 1069 jnoise = self.dataOut.getNoise()
511 1070 num_incoh = self.dataOut.nIncohInt
512 1071
513 1072 num_channel = jspectra.shape[0]
514 1073 num_prof = jspectra.shape[1]
515 1074 num_hei = jspectra.shape[2]
516 1075
517 1076 # hei_interf
518 1077 if hei_interf is None:
519 1078 count_hei = int(num_hei / 2)
520 1079 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
521 1080 hei_interf = numpy.asarray(hei_interf)[0]
522 1081 # nhei_interf
523 1082 if (nhei_interf == None):
524 1083 nhei_interf = 5
525 1084 if (nhei_interf < 1):
526 1085 nhei_interf = 1
527 1086 if (nhei_interf > count_hei):
528 1087 nhei_interf = count_hei
529 1088 if (offhei_interf == None):
530 1089 offhei_interf = 0
531 1090
532 1091 ind_hei = list(range(num_hei))
533 1092 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
534 1093 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
535 1094 mask_prof = numpy.asarray(list(range(num_prof)))
536 1095 num_mask_prof = mask_prof.size
537 1096 comp_mask_prof = [0, num_prof / 2]
538 1097
539 1098 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
540 1099 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
541 1100 jnoise = numpy.nan
542 1101 noise_exist = jnoise[0] < numpy.Inf
543 1102
544 1103 # Subrutina de Remocion de la Interferencia
545 1104 for ich in range(num_channel):
546 1105 # Se ordena los espectros segun su potencia (menor a mayor)
547 1106 power = jspectra[ich, mask_prof, :]
548 1107 power = power[:, hei_interf]
549 1108 power = power.sum(axis=0)
550 1109 psort = power.ravel().argsort()
551 1110
552 1111 # Se estima la interferencia promedio en los Espectros de Potencia empleando
553 1112 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
554 1113 offhei_interf, nhei_interf + offhei_interf))]]]
555 1114
556 1115 if noise_exist:
557 1116 # tmp_noise = jnoise[ich] / num_prof
558 1117 tmp_noise = jnoise[ich]
559 1118 junkspc_interf = junkspc_interf - tmp_noise
560 1119 #junkspc_interf[:,comp_mask_prof] = 0
561 1120
562 1121 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
563 1122 jspc_interf = jspc_interf.transpose()
564 1123 # Calculando el espectro de interferencia promedio
565 1124 noiseid = numpy.where(
566 1125 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
567 1126 noiseid = noiseid[0]
568 1127 cnoiseid = noiseid.size
569 1128 interfid = numpy.where(
570 1129 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
571 1130 interfid = interfid[0]
572 1131 cinterfid = interfid.size
573 1132
574 1133 if (cnoiseid > 0):
575 1134 jspc_interf[noiseid] = 0
576 1135
577 1136 # Expandiendo los perfiles a limpiar
578 1137 if (cinterfid > 0):
579 1138 new_interfid = (
580 1139 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
581 1140 new_interfid = numpy.asarray(new_interfid)
582 1141 new_interfid = {x for x in new_interfid}
583 1142 new_interfid = numpy.array(list(new_interfid))
584 1143 new_cinterfid = new_interfid.size
585 1144 else:
586 1145 new_cinterfid = 0
587 1146
588 1147 for ip in range(new_cinterfid):
589 1148 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
590 1149 jspc_interf[new_interfid[ip]
591 1150 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
592 1151
593 1152 jspectra[ich, :, ind_hei] = jspectra[ich, :,
594 1153 ind_hei] - jspc_interf # Corregir indices
595 1154
596 1155 # Removiendo la interferencia del punto de mayor interferencia
597 1156 ListAux = jspc_interf[mask_prof].tolist()
598 1157 maxid = ListAux.index(max(ListAux))
599 1158
600 1159 if cinterfid > 0:
601 1160 for ip in range(cinterfid * (interf == 2) - 1):
602 1161 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
603 1162 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
604 1163 cind = len(ind)
605 1164
606 1165 if (cind > 0):
607 1166 jspectra[ich, interfid[ip], ind] = tmp_noise * \
608 1167 (1 + (numpy.random.uniform(cind) - 0.5) /
609 1168 numpy.sqrt(num_incoh))
610 1169
611 1170 ind = numpy.array([-2, -1, 1, 2])
612 1171 xx = numpy.zeros([4, 4])
613 1172
614 1173 for id1 in range(4):
615 1174 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
616 1175
617 1176 xx_inv = numpy.linalg.inv(xx)
618 1177 xx = xx_inv[:, 0]
619 1178 ind = (ind + maxid + num_mask_prof) % num_mask_prof
620 1179 yy = jspectra[ich, mask_prof[ind], :]
621 1180 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
622 1181 yy.transpose(), xx)
623 1182
624 1183 indAux = (jspectra[ich, :, :] < tmp_noise *
625 1184 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
626 1185 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
627 1186 (1 - 1 / numpy.sqrt(num_incoh))
628 1187
629 1188 # Remocion de Interferencia en el Cross Spectra
630 1189 if jcspectra is None:
631 1190 return jspectra, jcspectra
632 1191 num_pairs = int(jcspectra.size / (num_prof * num_hei))
633 1192 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
634 1193
635 1194 for ip in range(num_pairs):
636 1195
637 1196 #-------------------------------------------
638 1197
639 1198 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
640 1199 cspower = cspower[:, hei_interf]
641 1200 cspower = cspower.sum(axis=0)
642 1201
643 1202 cspsort = cspower.ravel().argsort()
644 1203 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
645 1204 offhei_interf, nhei_interf + offhei_interf))]]]
646 1205 junkcspc_interf = junkcspc_interf.transpose()
647 1206 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
648 1207
649 1208 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
650 1209
651 1210 median_real = int(numpy.median(numpy.real(
652 1211 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
653 1212 median_imag = int(numpy.median(numpy.imag(
654 1213 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
655 1214 comp_mask_prof = [int(e) for e in comp_mask_prof]
656 1215 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
657 1216 median_real, median_imag)
658 1217
659 1218 for iprof in range(num_prof):
660 1219 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
661 1220 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
662 1221
663 1222 # Removiendo la Interferencia
664 1223 jcspectra[ip, :, ind_hei] = jcspectra[ip,
665 1224 :, ind_hei] - jcspc_interf
666 1225
667 1226 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
668 1227 maxid = ListAux.index(max(ListAux))
669 1228
670 1229 ind = numpy.array([-2, -1, 1, 2])
671 1230 xx = numpy.zeros([4, 4])
672 1231
673 1232 for id1 in range(4):
674 1233 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
675 1234
676 1235 xx_inv = numpy.linalg.inv(xx)
677 1236 xx = xx_inv[:, 0]
678 1237
679 1238 ind = (ind + maxid + num_mask_prof) % num_mask_prof
680 1239 yy = jcspectra[ip, mask_prof[ind], :]
681 1240 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
682 1241
683 1242 # Guardar Resultados
684 1243 self.dataOut.data_spc = jspectra
685 1244 self.dataOut.data_cspc = jcspectra
686 1245
687 1246 return 1
688 1247
689 1248 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
690 1249
691 1250 self.dataOut = dataOut
692 1251
693 1252 if mode == 1:
694 1253 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
695 1254 elif mode == 2:
696 1255 self.removeInterference2()
697 1256
698 1257 return self.dataOut
699 1258
700 1259
701 1260 class IncohInt(Operation):
702 1261
703 1262 __profIndex = 0
704 1263 __withOverapping = False
705 1264
706 1265 __byTime = False
707 1266 __initime = None
708 1267 __lastdatatime = None
709 1268 __integrationtime = None
710 1269
711 1270 __buffer_spc = None
712 1271 __buffer_cspc = None
713 1272 __buffer_dc = None
714 1273
715 1274 __dataReady = False
716 1275
717 1276 __timeInterval = None
718 1277
719 1278 n = None
720 1279
721 1280 def __init__(self):
722 1281
723 1282 Operation.__init__(self)
724 1283
725 1284 def setup(self, n=None, timeInterval=None, overlapping=False):
726 1285 """
727 1286 Set the parameters of the integration class.
728 1287
729 1288 Inputs:
730 1289
731 1290 n : Number of coherent integrations
732 1291 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
733 1292 overlapping :
734 1293
735 1294 """
736 1295
737 1296 self.__initime = None
738 1297 self.__lastdatatime = 0
739 1298
740 1299 self.__buffer_spc = 0
741 1300 self.__buffer_cspc = 0
742 1301 self.__buffer_dc = 0
743 1302
744 1303 self.__profIndex = 0
745 1304 self.__dataReady = False
746 1305 self.__byTime = False
747 1306
748 1307 if n is None and timeInterval is None:
749 1308 raise ValueError("n or timeInterval should be specified ...")
750 1309
751 1310 if n is not None:
752 1311 self.n = int(n)
753 1312 else:
754 1313
755 1314 self.__integrationtime = int(timeInterval)
756 1315 self.n = None
757 1316 self.__byTime = True
758 1317
759 1318 def putData(self, data_spc, data_cspc, data_dc):
760 1319 """
761 1320 Add a profile to the __buffer_spc and increase in one the __profileIndex
762 1321
763 1322 """
764 1323
765 1324 self.__buffer_spc += data_spc
766 1325
767 1326 if data_cspc is None:
768 1327 self.__buffer_cspc = None
769 1328 else:
770 1329 self.__buffer_cspc += data_cspc
771 1330
772 1331 if data_dc is None:
773 1332 self.__buffer_dc = None
774 1333 else:
775 1334 self.__buffer_dc += data_dc
776 1335
777 1336 self.__profIndex += 1
778 1337
779 1338 return
780 1339
781 1340 def pushData(self):
782 1341 """
783 1342 Return the sum of the last profiles and the profiles used in the sum.
784 1343
785 1344 Affected:
786 1345
787 1346 self.__profileIndex
788 1347
789 1348 """
790 1349
791 1350 data_spc = self.__buffer_spc
792 1351 data_cspc = self.__buffer_cspc
793 1352 data_dc = self.__buffer_dc
794 1353 n = self.__profIndex
795 1354
796 1355 self.__buffer_spc = 0
797 1356 self.__buffer_cspc = 0
798 1357 self.__buffer_dc = 0
799 1358 self.__profIndex = 0
800 1359
801 1360 return data_spc, data_cspc, data_dc, n
802 1361
803 1362 def byProfiles(self, *args):
804 1363
805 1364 self.__dataReady = False
806 1365 avgdata_spc = None
807 1366 avgdata_cspc = None
808 1367 avgdata_dc = None
809 1368
810 1369 self.putData(*args)
811 1370
812 1371 if self.__profIndex == self.n:
813 1372
814 1373 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
815 1374 self.n = n
816 1375 self.__dataReady = True
817 1376
818 1377 return avgdata_spc, avgdata_cspc, avgdata_dc
819 1378
820 1379 def byTime(self, datatime, *args):
821 1380
822 1381 self.__dataReady = False
823 1382 avgdata_spc = None
824 1383 avgdata_cspc = None
825 1384 avgdata_dc = None
826 1385
827 1386 self.putData(*args)
828 1387
829 1388 if (datatime - self.__initime) >= self.__integrationtime:
830 1389 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
831 1390 self.n = n
832 1391 self.__dataReady = True
833 1392
834 1393 return avgdata_spc, avgdata_cspc, avgdata_dc
835 1394
836 1395 def integrate(self, datatime, *args):
837 1396
838 1397 if self.__profIndex == 0:
839 1398 self.__initime = datatime
840 1399
841 1400 if self.__byTime:
842 1401 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
843 1402 datatime, *args)
844 1403 else:
845 1404 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
846 1405
847 1406 if not self.__dataReady:
848 1407 return None, None, None, None
849 1408
850 1409 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
851 1410
852 1411 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
853 1412 if n == 1:
854 1413 return dataOut
855 1414
856 1415 dataOut.flagNoData = True
857 1416
858 1417 if not self.isConfig:
859 1418 self.setup(n, timeInterval, overlapping)
860 1419 self.isConfig = True
861 1420
862 1421 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
863 1422 dataOut.data_spc,
864 1423 dataOut.data_cspc,
865 1424 dataOut.data_dc)
866 1425
867 1426 if self.__dataReady:
868 1427
869 1428 dataOut.data_spc = avgdata_spc
870 1429 dataOut.data_cspc = avgdata_cspc
871 1430 dataOut.data_dc = avgdata_dc
872 1431 dataOut.nIncohInt *= self.n
873 1432 dataOut.utctime = avgdatatime
874 1433 dataOut.flagNoData = False
875 1434
876 1435 return dataOut
877 1436
878 1437 class dopplerFlip(Operation):
879 1438
880 1439 def run(self, dataOut):
881 1440 # arreglo 1: (num_chan, num_profiles, num_heights)
882 1441 self.dataOut = dataOut
883 1442 # JULIA-oblicua, indice 2
884 1443 # arreglo 2: (num_profiles, num_heights)
885 1444 jspectra = self.dataOut.data_spc[2]
886 1445 jspectra_tmp = numpy.zeros(jspectra.shape)
887 1446 num_profiles = jspectra.shape[0]
888 1447 freq_dc = int(num_profiles / 2)
889 1448 # Flip con for
890 1449 for j in range(num_profiles):
891 1450 jspectra_tmp[num_profiles-j-1]= jspectra[j]
892 1451 # Intercambio perfil de DC con perfil inmediato anterior
893 1452 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
894 1453 jspectra_tmp[freq_dc]= jspectra[freq_dc]
895 1454 # canal modificado es re-escrito en el arreglo de canales
896 1455 self.dataOut.data_spc[2] = jspectra_tmp
897 1456
898 1457 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now