##// END OF EJS Templates
Multiprocessing for BLTR (sswma) data
Juan C. Espinoza -
r1181:eb60268405c6
parent child
Show More
@@ -1,369 +1,367
1 '''
1 '''
2 Created on Nov 9, 2016
2 Created on Nov 9, 2016
3
3
4 @author: roj- LouVD
4 @author: roj- LouVD
5 '''
5 '''
6
6
7
7
8 import os
8 import os
9 import sys
9 import sys
10 import time
10 import time
11 import glob
11 import glob
12 import datetime
12 import datetime
13
13
14 import numpy
14 import numpy
15
15
16 from schainpy.model.proc.jroproc_base import ProcessingUnit
16 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator
17 from schainpy.model.data.jrodata import Parameters
17 from schainpy.model.data.jrodata import Parameters
18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
18 from schainpy.model.io.jroIO_base import JRODataReader, isNumber
19 from schainpy.utils import log
19 from schainpy.utils import log
20
20
21 FILE_HEADER_STRUCTURE = numpy.dtype([
21 FILE_HEADER_STRUCTURE = numpy.dtype([
22 ('FMN', '<u4'),
22 ('FMN', '<u4'),
23 ('nrec', '<u4'),
23 ('nrec', '<u4'),
24 ('fr_offset', '<u4'),
24 ('fr_offset', '<u4'),
25 ('id', '<u4'),
25 ('id', '<u4'),
26 ('site', 'u1', (32,))
26 ('site', 'u1', (32,))
27 ])
27 ])
28
28
29 REC_HEADER_STRUCTURE = numpy.dtype([
29 REC_HEADER_STRUCTURE = numpy.dtype([
30 ('rmn', '<u4'),
30 ('rmn', '<u4'),
31 ('rcounter', '<u4'),
31 ('rcounter', '<u4'),
32 ('nr_offset', '<u4'),
32 ('nr_offset', '<u4'),
33 ('tr_offset', '<u4'),
33 ('tr_offset', '<u4'),
34 ('time', '<u4'),
34 ('time', '<u4'),
35 ('time_msec', '<u4'),
35 ('time_msec', '<u4'),
36 ('tag', 'u1', (32,)),
36 ('tag', 'u1', (32,)),
37 ('comments', 'u1', (32,)),
37 ('comments', 'u1', (32,)),
38 ('lat', '<f4'),
38 ('lat', '<f4'),
39 ('lon', '<f4'),
39 ('lon', '<f4'),
40 ('gps_status', '<u4'),
40 ('gps_status', '<u4'),
41 ('freq', '<u4'),
41 ('freq', '<u4'),
42 ('freq0', '<u4'),
42 ('freq0', '<u4'),
43 ('nchan', '<u4'),
43 ('nchan', '<u4'),
44 ('delta_r', '<u4'),
44 ('delta_r', '<u4'),
45 ('nranges', '<u4'),
45 ('nranges', '<u4'),
46 ('r0', '<u4'),
46 ('r0', '<u4'),
47 ('prf', '<u4'),
47 ('prf', '<u4'),
48 ('ncoh', '<u4'),
48 ('ncoh', '<u4'),
49 ('npoints', '<u4'),
49 ('npoints', '<u4'),
50 ('polarization', '<i4'),
50 ('polarization', '<i4'),
51 ('rx_filter', '<u4'),
51 ('rx_filter', '<u4'),
52 ('nmodes', '<u4'),
52 ('nmodes', '<u4'),
53 ('dmode_index', '<u4'),
53 ('dmode_index', '<u4'),
54 ('dmode_rngcorr', '<u4'),
54 ('dmode_rngcorr', '<u4'),
55 ('nrxs', '<u4'),
55 ('nrxs', '<u4'),
56 ('acf_length', '<u4'),
56 ('acf_length', '<u4'),
57 ('acf_lags', '<u4'),
57 ('acf_lags', '<u4'),
58 ('sea_to_atmos', '<f4'),
58 ('sea_to_atmos', '<f4'),
59 ('sea_notch', '<u4'),
59 ('sea_notch', '<u4'),
60 ('lh_sea', '<u4'),
60 ('lh_sea', '<u4'),
61 ('hh_sea', '<u4'),
61 ('hh_sea', '<u4'),
62 ('nbins_sea', '<u4'),
62 ('nbins_sea', '<u4'),
63 ('min_snr', '<f4'),
63 ('min_snr', '<f4'),
64 ('min_cc', '<f4'),
64 ('min_cc', '<f4'),
65 ('max_time_diff', '<f4')
65 ('max_time_diff', '<f4')
66 ])
66 ])
67
67
68 DATA_STRUCTURE = numpy.dtype([
68 DATA_STRUCTURE = numpy.dtype([
69 ('range', '<u4'),
69 ('range', '<u4'),
70 ('status', '<u4'),
70 ('status', '<u4'),
71 ('zonal', '<f4'),
71 ('zonal', '<f4'),
72 ('meridional', '<f4'),
72 ('meridional', '<f4'),
73 ('vertical', '<f4'),
73 ('vertical', '<f4'),
74 ('zonal_a', '<f4'),
74 ('zonal_a', '<f4'),
75 ('meridional_a', '<f4'),
75 ('meridional_a', '<f4'),
76 ('corrected_fading', '<f4'), # seconds
76 ('corrected_fading', '<f4'), # seconds
77 ('uncorrected_fading', '<f4'), # seconds
77 ('uncorrected_fading', '<f4'), # seconds
78 ('time_diff', '<f4'),
78 ('time_diff', '<f4'),
79 ('major_axis', '<f4'),
79 ('major_axis', '<f4'),
80 ('axial_ratio', '<f4'),
80 ('axial_ratio', '<f4'),
81 ('orientation', '<f4'),
81 ('orientation', '<f4'),
82 ('sea_power', '<u4'),
82 ('sea_power', '<u4'),
83 ('sea_algorithm', '<u4')
83 ('sea_algorithm', '<u4')
84 ])
84 ])
85
85
86
86 @MPDecorator
87 class BLTRParamReader(JRODataReader, ProcessingUnit):
87 class BLTRParamReader(JRODataReader, ProcessingUnit):
88 '''
88 '''
89 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
89 Boundary Layer and Tropospheric Radar (BLTR) reader, Wind velocities and SNR from *.sswma files
90 '''
90 '''
91
91
92 ext = '.sswma'
92 ext = '.sswma'
93
93
94 def __init__(self, **kwargs):
94 def __init__(self, **kwargs):
95
95
96 ProcessingUnit.__init__(self, **kwargs)
96 ProcessingUnit.__init__(self, **kwargs)
97
97
98 self.dataOut = Parameters()
98 self.dataOut = Parameters()
99 self.counter_records = 0
99 self.counter_records = 0
100 self.flagNoMoreFiles = 0
100 self.flagNoMoreFiles = 0
101 self.isConfig = False
101 self.isConfig = False
102 self.filename = None
102 self.filename = None
103
103
104 def setup(self,
104 def setup(self,
105 path=None,
105 path=None,
106 startDate=None,
106 startDate=None,
107 endDate=None,
107 endDate=None,
108 ext=None,
108 ext=None,
109 startTime=datetime.time(0, 0, 0),
109 startTime=datetime.time(0, 0, 0),
110 endTime=datetime.time(23, 59, 59),
110 endTime=datetime.time(23, 59, 59),
111 timezone=0,
111 timezone=0,
112 status_value=0,
112 status_value=0,
113 **kwargs):
113 **kwargs):
114
114
115 self.path = path
115 self.path = path
116 self.startDate = startDate
116 self.startDate = startDate
117 self.endDate = endDate
117 self.endDate = endDate
118 self.startTime = startTime
118 self.startTime = startTime
119 self.endTime = endTime
119 self.endTime = endTime
120 self.status_value = status_value
120 self.status_value = status_value
121 self.datatime = datetime.datetime(1900,1,1)
121 self.datatime = datetime.datetime(1900,1,1)
122
122
123 if self.path is None:
123 if self.path is None:
124 raise ValueError("The path is not valid")
124 raise ValueError("The path is not valid")
125
125
126 if ext is None:
126 if ext is None:
127 ext = self.ext
127 ext = self.ext
128
128
129 self.search_files(self.path, startDate, endDate, ext)
129 self.search_files(self.path, startDate, endDate, ext)
130 self.timezone = timezone
130 self.timezone = timezone
131 self.fileIndex = 0
131 self.fileIndex = 0
132
132
133 if not self.fileList:
133 if not self.fileList:
134 raise Warning("There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' " % (
134 raise Warning("There is no files matching these date in the folder: %s. \n Check 'startDate' and 'endDate' " % (
135 path))
135 path))
136
136
137 self.setNextFile()
137 self.setNextFile()
138
138
139 def search_files(self, path, startDate, endDate, ext):
139 def search_files(self, path, startDate, endDate, ext):
140 '''
140 '''
141 Searching for BLTR rawdata file in path
141 Searching for BLTR rawdata file in path
142 Creating a list of file to proces included in [startDate,endDate]
142 Creating a list of file to proces included in [startDate,endDate]
143
143
144 Input:
144 Input:
145 path - Path to find BLTR rawdata files
145 path - Path to find BLTR rawdata files
146 startDate - Select file from this date
146 startDate - Select file from this date
147 enDate - Select file until this date
147 enDate - Select file until this date
148 ext - Extension of the file to read
148 ext - Extension of the file to read
149 '''
149 '''
150
150
151 log.success('Searching files in {} '.format(path), 'BLTRParamReader')
151 log.success('Searching files in {} '.format(path), 'BLTRParamReader')
152 foldercounter = 0
152 foldercounter = 0
153 fileList0 = glob.glob1(path, "*%s" % ext)
153 fileList0 = glob.glob1(path, "*%s" % ext)
154 fileList0.sort()
154 fileList0.sort()
155
155
156 self.fileList = []
156 self.fileList = []
157 self.dateFileList = []
157 self.dateFileList = []
158
158
159 for thisFile in fileList0:
159 for thisFile in fileList0:
160 year = thisFile[-14:-10]
160 year = thisFile[-14:-10]
161 if not isNumber(year):
161 if not isNumber(year):
162 continue
162 continue
163
163
164 month = thisFile[-10:-8]
164 month = thisFile[-10:-8]
165 if not isNumber(month):
165 if not isNumber(month):
166 continue
166 continue
167
167
168 day = thisFile[-8:-6]
168 day = thisFile[-8:-6]
169 if not isNumber(day):
169 if not isNumber(day):
170 continue
170 continue
171
171
172 year, month, day = int(year), int(month), int(day)
172 year, month, day = int(year), int(month), int(day)
173 dateFile = datetime.date(year, month, day)
173 dateFile = datetime.date(year, month, day)
174
174
175 if (startDate > dateFile) or (endDate < dateFile):
175 if (startDate > dateFile) or (endDate < dateFile):
176 continue
176 continue
177
177
178 self.fileList.append(thisFile)
178 self.fileList.append(thisFile)
179 self.dateFileList.append(dateFile)
179 self.dateFileList.append(dateFile)
180
180
181 return
181 return
182
182
183 def setNextFile(self):
183 def setNextFile(self):
184
184
185 file_id = self.fileIndex
185 file_id = self.fileIndex
186
186
187 if file_id == len(self.fileList):
187 if file_id == len(self.fileList):
188 log.success('No more files in the folder', 'BLTRParamReader')
189 self.flagNoMoreFiles = 1
188 self.flagNoMoreFiles = 1
190 return 0
189 return 0
191
190
192 log.success('Opening {}'.format(self.fileList[file_id]), 'BLTRParamReader')
191 log.success('Opening {}'.format(self.fileList[file_id]), 'BLTRParamReader')
193 filename = os.path.join(self.path, self.fileList[file_id])
192 filename = os.path.join(self.path, self.fileList[file_id])
194
193
195 dirname, name = os.path.split(filename)
194 dirname, name = os.path.split(filename)
196 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
195 # 'peru2' ---> Piura - 'peru1' ---> Huancayo or Porcuya
197 self.siteFile = name.split('.')[0]
196 self.siteFile = name.split('.')[0]
198 if self.filename is not None:
197 if self.filename is not None:
199 self.fp.close()
198 self.fp.close()
200 self.filename = filename
199 self.filename = filename
201 self.fp = open(self.filename, 'rb')
200 self.fp = open(self.filename, 'rb')
202 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
201 self.header_file = numpy.fromfile(self.fp, FILE_HEADER_STRUCTURE, 1)
203 self.nrecords = self.header_file['nrec'][0]
202 self.nrecords = self.header_file['nrec'][0]
204 self.sizeOfFile = os.path.getsize(self.filename)
203 self.sizeOfFile = os.path.getsize(self.filename)
205 self.counter_records = 0
204 self.counter_records = 0
206 self.flagIsNewFile = 0
205 self.flagIsNewFile = 0
207 self.fileIndex += 1
206 self.fileIndex += 1
208
207
209 return 1
208 return 1
210
209
211 def readNextBlock(self):
210 def readNextBlock(self):
212
211
213 while True:
212 while True:
214 if self.counter_records == self.nrecords:
213 if self.counter_records == self.nrecords:
215 self.flagIsNewFile = 1
214 self.flagIsNewFile = 1
216 if not self.setNextFile():
215 if not self.setNextFile():
217 return 0
216 return 0
218
217
219 self.readBlock()
218 self.readBlock()
220
219
221 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
220 if (self.datatime < datetime.datetime.combine(self.startDate, self.startTime)) or \
222 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
221 (self.datatime > datetime.datetime.combine(self.endDate, self.endTime)):
223 log.warning(
222 log.warning(
224 'Reading Record No. {}/{} -> {} [Skipping]'.format(
223 'Reading Record No. {}/{} -> {} [Skipping]'.format(
225 self.counter_records,
224 self.counter_records,
226 self.nrecords,
225 self.nrecords,
227 self.datatime.ctime()),
226 self.datatime.ctime()),
228 'BLTRParamReader')
227 'BLTRParamReader')
229 continue
228 continue
230 break
229 break
231
230
232 log.log('Reading Record No. {}/{} -> {}'.format(
231 log.log('Reading Record No. {}/{} -> {}'.format(
233 self.counter_records,
232 self.counter_records,
234 self.nrecords,
233 self.nrecords,
235 self.datatime.ctime()), 'BLTRParamReader')
234 self.datatime.ctime()), 'BLTRParamReader')
236
235
237 return 1
236 return 1
238
237
239 def readBlock(self):
238 def readBlock(self):
240
239
241 pointer = self.fp.tell()
240 pointer = self.fp.tell()
242 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
241 header_rec = numpy.fromfile(self.fp, REC_HEADER_STRUCTURE, 1)
243 self.nchannels = header_rec['nchan'][0] / 2
242 self.nchannels = int(header_rec['nchan'][0] / 2)
244 self.kchan = header_rec['nrxs'][0]
243 self.kchan = header_rec['nrxs'][0]
245 self.nmodes = header_rec['nmodes'][0]
244 self.nmodes = header_rec['nmodes'][0]
246 self.nranges = header_rec['nranges'][0]
245 self.nranges = header_rec['nranges'][0]
247 self.fp.seek(pointer)
246 self.fp.seek(pointer)
248 self.height = numpy.empty((self.nmodes, self.nranges))
247 self.height = numpy.empty((self.nmodes, self.nranges))
249 self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
248 self.snr = numpy.empty((self.nmodes, self.nchannels, self.nranges))
250 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
249 self.buffer = numpy.empty((self.nmodes, 3, self.nranges))
251 self.flagDiscontinuousBlock = 0
250 self.flagDiscontinuousBlock = 0
252
251
253 for mode in range(self.nmodes):
252 for mode in range(self.nmodes):
254 self.readHeader()
253 self.readHeader()
255 data = self.readData()
254 data = self.readData()
256 self.height[mode] = (data[0] - self.correction) / 1000.
255 self.height[mode] = (data[0] - self.correction) / 1000.
257 self.buffer[mode] = data[1]
256 self.buffer[mode] = data[1]
258 self.snr[mode] = data[2]
257 self.snr[mode] = data[2]
259
258
260 self.counter_records = self.counter_records + self.nmodes
259 self.counter_records = self.counter_records + self.nmodes
261
260
262 return
261 return
263
262
264 def readHeader(self):
263 def readHeader(self):
265 '''
264 '''
266 RecordHeader of BLTR rawdata file
265 RecordHeader of BLTR rawdata file
267 '''
266 '''
268
267
269 header_structure = numpy.dtype(
268 header_structure = numpy.dtype(
270 REC_HEADER_STRUCTURE.descr + [
269 REC_HEADER_STRUCTURE.descr + [
271 ('antenna_coord', 'f4', (2, self.nchannels)),
270 ('antenna_coord', 'f4', (2, self.nchannels)),
272 ('rx_gains', 'u4', (self.nchannels,)),
271 ('rx_gains', 'u4', (self.nchannels,)),
273 ('rx_analysis', 'u4', (self.nchannels,))
272 ('rx_analysis', 'u4', (self.nchannels,))
274 ]
273 ]
275 )
274 )
276
275
277 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
276 self.header_rec = numpy.fromfile(self.fp, header_structure, 1)
278 self.lat = self.header_rec['lat'][0]
277 self.lat = self.header_rec['lat'][0]
279 self.lon = self.header_rec['lon'][0]
278 self.lon = self.header_rec['lon'][0]
280 self.delta = self.header_rec['delta_r'][0]
279 self.delta = self.header_rec['delta_r'][0]
281 self.correction = self.header_rec['dmode_rngcorr'][0]
280 self.correction = self.header_rec['dmode_rngcorr'][0]
282 self.imode = self.header_rec['dmode_index'][0]
281 self.imode = self.header_rec['dmode_index'][0]
283 self.antenna = self.header_rec['antenna_coord']
282 self.antenna = self.header_rec['antenna_coord']
284 self.rx_gains = self.header_rec['rx_gains']
283 self.rx_gains = self.header_rec['rx_gains']
285 self.time = self.header_rec['time'][0]
284 self.time = self.header_rec['time'][0]
286 dt = datetime.datetime.utcfromtimestamp(self.time)
285 dt = datetime.datetime.utcfromtimestamp(self.time)
287 if dt.date()>self.datatime.date():
286 if dt.date()>self.datatime.date():
288 self.flagDiscontinuousBlock = 1
287 self.flagDiscontinuousBlock = 1
289 self.datatime = dt
288 self.datatime = dt
290
289
291 def readData(self):
290 def readData(self):
292 '''
291 '''
293 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
292 Reading and filtering data block record of BLTR rawdata file, filtering is according to status_value.
294
293
295 Input:
294 Input:
296 status_value - Array data is set to NAN for values that are not equal to status_value
295 status_value - Array data is set to NAN for values that are not equal to status_value
297
296
298 '''
297 '''
299
298
300 data_structure = numpy.dtype(
299 data_structure = numpy.dtype(
301 DATA_STRUCTURE.descr + [
300 DATA_STRUCTURE.descr + [
302 ('rx_saturation', 'u4', (self.nchannels,)),
301 ('rx_saturation', 'u4', (self.nchannels,)),
303 ('chan_offset', 'u4', (2 * self.nchannels,)),
302 ('chan_offset', 'u4', (2 * self.nchannels,)),
304 ('rx_amp', 'u4', (self.nchannels,)),
303 ('rx_amp', 'u4', (self.nchannels,)),
305 ('rx_snr', 'f4', (self.nchannels,)),
304 ('rx_snr', 'f4', (self.nchannels,)),
306 ('cross_snr', 'f4', (self.kchan,)),
305 ('cross_snr', 'f4', (self.kchan,)),
307 ('sea_power_relative', 'f4', (self.kchan,))]
306 ('sea_power_relative', 'f4', (self.kchan,))]
308 )
307 )
309
308
310 data = numpy.fromfile(self.fp, data_structure, self.nranges)
309 data = numpy.fromfile(self.fp, data_structure, self.nranges)
311
310
312 height = data['range']
311 height = data['range']
313 winds = numpy.array(
312 winds = numpy.array(
314 (data['zonal'], data['meridional'], data['vertical']))
313 (data['zonal'], data['meridional'], data['vertical']))
315 snr = data['rx_snr'].T
314 snr = data['rx_snr'].T
316
315
317 winds[numpy.where(winds == -9999.)] = numpy.nan
316 winds[numpy.where(winds == -9999.)] = numpy.nan
318 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
317 winds[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
319 snr[numpy.where(snr == -9999.)] = numpy.nan
318 snr[numpy.where(snr == -9999.)] = numpy.nan
320 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
319 snr[:, numpy.where(data['status'] != self.status_value)] = numpy.nan
321 snr = numpy.power(10, snr / 10)
320 snr = numpy.power(10, snr / 10)
322
321
323 return height, winds, snr
322 return height, winds, snr
324
323
325 def set_output(self):
324 def set_output(self):
326 '''
325 '''
327 Storing data from databuffer to dataOut object
326 Storing data from databuffer to dataOut object
328 '''
327 '''
329
328
330 self.dataOut.data_SNR = self.snr
329 self.dataOut.data_SNR = self.snr
331 self.dataOut.height = self.height
330 self.dataOut.height = self.height
332 self.dataOut.data = self.buffer
331 self.dataOut.data = self.buffer
333 self.dataOut.utctimeInit = self.time
332 self.dataOut.utctimeInit = self.time
334 self.dataOut.utctime = self.dataOut.utctimeInit
333 self.dataOut.utctime = self.dataOut.utctimeInit
335 self.dataOut.useLocalTime = False
334 self.dataOut.useLocalTime = False
336 self.dataOut.paramInterval = 157
335 self.dataOut.paramInterval = 157
337 self.dataOut.timezone = self.timezone
336 self.dataOut.timezone = self.timezone
338 self.dataOut.site = self.siteFile
337 self.dataOut.site = self.siteFile
339 self.dataOut.nrecords = self.nrecords / self.nmodes
338 self.dataOut.nrecords = self.nrecords / self.nmodes
340 self.dataOut.sizeOfFile = self.sizeOfFile
339 self.dataOut.sizeOfFile = self.sizeOfFile
341 self.dataOut.lat = self.lat
340 self.dataOut.lat = self.lat
342 self.dataOut.lon = self.lon
341 self.dataOut.lon = self.lon
343 self.dataOut.channelList = list(range(self.nchannels))
342 self.dataOut.channelList = list(range(self.nchannels))
344 self.dataOut.kchan = self.kchan
343 self.dataOut.kchan = self.kchan
345 self.dataOut.delta = self.delta
344 self.dataOut.delta = self.delta
346 self.dataOut.correction = self.correction
345 self.dataOut.correction = self.correction
347 self.dataOut.nmodes = self.nmodes
346 self.dataOut.nmodes = self.nmodes
348 self.dataOut.imode = self.imode
347 self.dataOut.imode = self.imode
349 self.dataOut.antenna = self.antenna
348 self.dataOut.antenna = self.antenna
350 self.dataOut.rx_gains = self.rx_gains
349 self.dataOut.rx_gains = self.rx_gains
351 self.dataOut.flagNoData = False
350 self.dataOut.flagNoData = False
352 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
351 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
353
352
354 def getData(self):
353 def getData(self):
355 '''
354 '''
356 Storing data from databuffer to dataOut object
355 Storing data from databuffer to dataOut object
357 '''
356 '''
358 if self.flagNoMoreFiles:
357 if self.flagNoMoreFiles:
359 self.dataOut.flagNoData = True
358 self.dataOut.flagNoData = True
360 log.success('No file left to process', 'BLTRParamReader')
359 self.dataOut.error = (1, 'No More files to read')
361 return 0
362
360
363 if not self.readNextBlock():
361 if not self.readNextBlock():
364 self.dataOut.flagNoData = True
362 self.dataOut.flagNoData = True
365 return 0
363 return 0
366
364
367 self.set_output()
365 self.set_output()
368
366
369 return 1 No newline at end of file
367 return 1
@@ -1,402 +1,401
1 '''
1 '''
2 Created on Oct 24, 2016
2 Created on Oct 24, 2016
3
3
4 @author: roj- LouVD
4 @author: roj- LouVD
5 '''
5 '''
6
6
7 import numpy
7 import numpy
8 import copy
8 import copy
9 import datetime
9 import datetime
10 import time
10 import time
11 from time import gmtime
11 from time import gmtime
12
12
13 from numpy import transpose
13 from numpy import transpose
14
14
15 from .jroproc_base import ProcessingUnit, Operation
15 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
16 from schainpy.model.data.jrodata import Parameters
16 from schainpy.model.data.jrodata import Parameters
17
17
18
18 @MPDecorator
19 class BLTRParametersProc(ProcessingUnit):
19 class BLTRParametersProc(ProcessingUnit):
20 '''
20 '''
21 Processing unit for BLTR parameters data (winds)
21 Processing unit for BLTR parameters data (winds)
22
22
23 Inputs:
23 Inputs:
24 self.dataOut.nmodes - Number of operation modes
24 self.dataOut.nmodes - Number of operation modes
25 self.dataOut.nchannels - Number of channels
25 self.dataOut.nchannels - Number of channels
26 self.dataOut.nranges - Number of ranges
26 self.dataOut.nranges - Number of ranges
27
27
28 self.dataOut.data_SNR - SNR array
28 self.dataOut.data_SNR - SNR array
29 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
29 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
30 self.dataOut.height - Height array (km)
30 self.dataOut.height - Height array (km)
31 self.dataOut.time - Time array (seconds)
31 self.dataOut.time - Time array (seconds)
32
32
33 self.dataOut.fileIndex -Index of the file currently read
33 self.dataOut.fileIndex -Index of the file currently read
34 self.dataOut.lat - Latitude coordinate of BLTR location
34 self.dataOut.lat - Latitude coordinate of BLTR location
35
35
36 self.dataOut.doy - Experiment doy (number of the day in the current year)
36 self.dataOut.doy - Experiment doy (number of the day in the current year)
37 self.dataOut.month - Experiment month
37 self.dataOut.month - Experiment month
38 self.dataOut.day - Experiment day
38 self.dataOut.day - Experiment day
39 self.dataOut.year - Experiment year
39 self.dataOut.year - Experiment year
40 '''
40 '''
41
41
42 def __init__(self, **kwargs):
42 def __init__(self):
43 '''
43 '''
44 Inputs: None
44 Inputs: None
45 '''
45 '''
46 ProcessingUnit.__init__(self, **kwargs)
46 ProcessingUnit.__init__(self)
47 self.dataOut = Parameters()
47 self.dataOut = Parameters()
48 self.isConfig = False
49
48
50 def setup(self, mode):
49 def setup(self, mode):
51 '''
50 '''
52 '''
51 '''
53 self.dataOut.mode = mode
52 self.dataOut.mode = mode
54
53
55 def run(self, mode, snr_threshold=None):
54 def run(self, mode, snr_threshold=None):
56 '''
55 '''
57 Inputs:
56 Inputs:
58 mode = High resolution (0) or Low resolution (1) data
57 mode = High resolution (0) or Low resolution (1) data
59 snr_threshold = snr filter value
58 snr_threshold = snr filter value
60 '''
59 '''
61
60
62 if not self.isConfig:
61 if not self.isConfig:
63 self.setup(mode)
62 self.setup(mode)
64 self.isConfig = True
63 self.isConfig = True
65
64
66 if self.dataIn.type == 'Parameters':
65 if self.dataIn.type == 'Parameters':
67 self.dataOut.copy(self.dataIn)
66 self.dataOut.copy(self.dataIn)
68
67
69 self.dataOut.data_param = self.dataOut.data[mode]
68 self.dataOut.data_param = self.dataOut.data[mode]
70 self.dataOut.heightList = self.dataOut.height[0]
69 self.dataOut.heightList = self.dataOut.height[0]
71 self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
70 self.dataOut.data_SNR = self.dataOut.data_SNR[mode]
72
71
73 if snr_threshold is not None:
72 if snr_threshold is not None:
74 SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
73 SNRavg = numpy.average(self.dataOut.data_SNR, axis=0)
75 SNRavgdB = 10*numpy.log10(SNRavg)
74 SNRavgdB = 10*numpy.log10(SNRavg)
76 for i in range(3):
75 for i in range(3):
77 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
76 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
78
77
79 # TODO
78 # TODO
80 class OutliersFilter(Operation):
79 class OutliersFilter(Operation):
81
80
82 def __init__(self, **kwargs):
81 def __init__(self, **kwargs):
83 '''
82 '''
84 '''
83 '''
85 Operation.__init__(self, **kwargs)
84 Operation.__init__(self, **kwargs)
86
85
87 def run(self, svalue2, method, factor, filter, npoints=9):
86 def run(self, svalue2, method, factor, filter, npoints=9):
88 '''
87 '''
89 Inputs:
88 Inputs:
90 svalue - string to select array velocity
89 svalue - string to select array velocity
91 svalue2 - string to choose axis filtering
90 svalue2 - string to choose axis filtering
92 method - 0 for SMOOTH or 1 for MEDIAN
91 method - 0 for SMOOTH or 1 for MEDIAN
93 factor - number used to set threshold
92 factor - number used to set threshold
94 filter - 1 for data filtering using the standard deviation criteria else 0
93 filter - 1 for data filtering using the standard deviation criteria else 0
95 npoints - number of points for mask filter
94 npoints - number of points for mask filter
96 '''
95 '''
97
96
98 print(' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor))
97 print(' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor))
99
98
100
99
101 yaxis = self.dataOut.heightList
100 yaxis = self.dataOut.heightList
102 xaxis = numpy.array([[self.dataOut.utctime]])
101 xaxis = numpy.array([[self.dataOut.utctime]])
103
102
104 # Zonal
103 # Zonal
105 value_temp = self.dataOut.data_output[0]
104 value_temp = self.dataOut.data_output[0]
106
105
107 # Zonal
106 # Zonal
108 value_temp = self.dataOut.data_output[1]
107 value_temp = self.dataOut.data_output[1]
109
108
110 # Vertical
109 # Vertical
111 value_temp = numpy.transpose(self.dataOut.data_output[2])
110 value_temp = numpy.transpose(self.dataOut.data_output[2])
112
111
113 htemp = yaxis
112 htemp = yaxis
114 std = value_temp
113 std = value_temp
115 for h in range(len(htemp)):
114 for h in range(len(htemp)):
116 nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
115 nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
117 minvalid = npoints
116 minvalid = npoints
118
117
119 #only if valid values greater than the minimum required (10%)
118 #only if valid values greater than the minimum required (10%)
120 if nvalues_valid > minvalid:
119 if nvalues_valid > minvalid:
121
120
122 if method == 0:
121 if method == 0:
123 #SMOOTH
122 #SMOOTH
124 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
123 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
125
124
126
125
127 if method == 1:
126 if method == 1:
128 #MEDIAN
127 #MEDIAN
129 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
128 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
130
129
131 dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
130 dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
132
131
133 threshold = dw*factor
132 threshold = dw*factor
134 value_temp[numpy.where(w > threshold),h] = numpy.nan
133 value_temp[numpy.where(w > threshold),h] = numpy.nan
135 value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
134 value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
136
135
137
136
138 #At the end
137 #At the end
139 if svalue2 == 'inHeight':
138 if svalue2 == 'inHeight':
140 value_temp = numpy.transpose(value_temp)
139 value_temp = numpy.transpose(value_temp)
141 output_array[:,m] = value_temp
140 output_array[:,m] = value_temp
142
141
143 if svalue == 'zonal':
142 if svalue == 'zonal':
144 self.dataOut.data_output[0] = output_array
143 self.dataOut.data_output[0] = output_array
145
144
146 elif svalue == 'meridional':
145 elif svalue == 'meridional':
147 self.dataOut.data_output[1] = output_array
146 self.dataOut.data_output[1] = output_array
148
147
149 elif svalue == 'vertical':
148 elif svalue == 'vertical':
150 self.dataOut.data_output[2] = output_array
149 self.dataOut.data_output[2] = output_array
151
150
152 return self.dataOut.data_output
151 return self.dataOut.data_output
153
152
154
153
155 def Median(self,input,width):
154 def Median(self,input,width):
156 '''
155 '''
157 Inputs:
156 Inputs:
158 input - Velocity array
157 input - Velocity array
159 width - Number of points for mask filter
158 width - Number of points for mask filter
160
159
161 '''
160 '''
162
161
163 if numpy.mod(width,2) == 1:
162 if numpy.mod(width,2) == 1:
164 pc = int((width - 1) / 2)
163 pc = int((width - 1) / 2)
165 cont = 0
164 cont = 0
166 output = []
165 output = []
167
166
168 for i in range(len(input)):
167 for i in range(len(input)):
169 if i >= pc and i < len(input) - pc:
168 if i >= pc and i < len(input) - pc:
170 new2 = input[i-pc:i+pc+1]
169 new2 = input[i-pc:i+pc+1]
171 temp = numpy.where(numpy.isfinite(new2))
170 temp = numpy.where(numpy.isfinite(new2))
172 new = new2[temp]
171 new = new2[temp]
173 value = numpy.median(new)
172 value = numpy.median(new)
174 output.append(value)
173 output.append(value)
175
174
176 output = numpy.array(output)
175 output = numpy.array(output)
177 output = numpy.hstack((input[0:pc],output))
176 output = numpy.hstack((input[0:pc],output))
178 output = numpy.hstack((output,input[-pc:len(input)]))
177 output = numpy.hstack((output,input[-pc:len(input)]))
179
178
180 return output
179 return output
181
180
182 def Smooth(self,input,width,edge_truncate = None):
181 def Smooth(self,input,width,edge_truncate = None):
183 '''
182 '''
184 Inputs:
183 Inputs:
185 input - Velocity array
184 input - Velocity array
186 width - Number of points for mask filter
185 width - Number of points for mask filter
187 edge_truncate - 1 for truncate the convolution product else
186 edge_truncate - 1 for truncate the convolution product else
188
187
189 '''
188 '''
190
189
191 if numpy.mod(width,2) == 0:
190 if numpy.mod(width,2) == 0:
192 real_width = width + 1
191 real_width = width + 1
193 nzeros = width / 2
192 nzeros = width / 2
194 else:
193 else:
195 real_width = width
194 real_width = width
196 nzeros = (width - 1) / 2
195 nzeros = (width - 1) / 2
197
196
198 half_width = int(real_width)/2
197 half_width = int(real_width)/2
199 length = len(input)
198 length = len(input)
200
199
201 gate = numpy.ones(real_width,dtype='float')
200 gate = numpy.ones(real_width,dtype='float')
202 norm_of_gate = numpy.sum(gate)
201 norm_of_gate = numpy.sum(gate)
203
202
204 nan_process = 0
203 nan_process = 0
205 nan_id = numpy.where(numpy.isnan(input))
204 nan_id = numpy.where(numpy.isnan(input))
206 if len(nan_id[0]) > 0:
205 if len(nan_id[0]) > 0:
207 nan_process = 1
206 nan_process = 1
208 pb = numpy.zeros(len(input))
207 pb = numpy.zeros(len(input))
209 pb[nan_id] = 1.
208 pb[nan_id] = 1.
210 input[nan_id] = 0.
209 input[nan_id] = 0.
211
210
212 if edge_truncate == True:
211 if edge_truncate == True:
213 output = numpy.convolve(input/norm_of_gate,gate,mode='same')
212 output = numpy.convolve(input/norm_of_gate,gate,mode='same')
214 elif edge_truncate == False or edge_truncate == None:
213 elif edge_truncate == False or edge_truncate == None:
215 output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
214 output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
216 output = numpy.hstack((input[0:half_width],output))
215 output = numpy.hstack((input[0:half_width],output))
217 output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
216 output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
218
217
219 if nan_process:
218 if nan_process:
220 pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
219 pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
221 pb = numpy.hstack((numpy.zeros(half_width),pb))
220 pb = numpy.hstack((numpy.zeros(half_width),pb))
222 pb = numpy.hstack((pb,numpy.zeros(half_width)))
221 pb = numpy.hstack((pb,numpy.zeros(half_width)))
223 output[numpy.where(pb > 0.9999)] = numpy.nan
222 output[numpy.where(pb > 0.9999)] = numpy.nan
224 input[nan_id] = numpy.nan
223 input[nan_id] = numpy.nan
225 return output
224 return output
226
225
227 def Average(self,aver=0,nhaver=1):
226 def Average(self,aver=0,nhaver=1):
228 '''
227 '''
229 Inputs:
228 Inputs:
230 aver - Indicates the time period over which is averaged or consensus data
229 aver - Indicates the time period over which is averaged or consensus data
231 nhaver - Indicates the decimation factor in heights
230 nhaver - Indicates the decimation factor in heights
232
231
233 '''
232 '''
234 nhpoints = 48
233 nhpoints = 48
235
234
236 lat_piura = -5.17
235 lat_piura = -5.17
237 lat_huancayo = -12.04
236 lat_huancayo = -12.04
238 lat_porcuya = -5.8
237 lat_porcuya = -5.8
239
238
240 if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
239 if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
241 hcm = 3.
240 hcm = 3.
242 if self.dataOut.year == 2003 :
241 if self.dataOut.year == 2003 :
243 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
242 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
244 nhpoints = 12
243 nhpoints = 12
245
244
246 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
245 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
247 hcm = 3.
246 hcm = 3.
248 if self.dataOut.year == 2003 :
247 if self.dataOut.year == 2003 :
249 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
248 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
250 nhpoints = 12
249 nhpoints = 12
251
250
252
251
253 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
252 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
254 hcm = 5.#2
253 hcm = 5.#2
255
254
256 pdata = 0.2
255 pdata = 0.2
257 taver = [1,2,3,4,6,8,12,24]
256 taver = [1,2,3,4,6,8,12,24]
258 t0 = 0
257 t0 = 0
259 tf = 24
258 tf = 24
260 ntime =(tf-t0)/taver[aver]
259 ntime =(tf-t0)/taver[aver]
261 ti = numpy.arange(ntime)
260 ti = numpy.arange(ntime)
262 tf = numpy.arange(ntime) + taver[aver]
261 tf = numpy.arange(ntime) + taver[aver]
263
262
264
263
265 old_height = self.dataOut.heightList
264 old_height = self.dataOut.heightList
266
265
267 if nhaver > 1:
266 if nhaver > 1:
268 num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
267 num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
269 deltha = 0.05*nhaver
268 deltha = 0.05*nhaver
270 minhvalid = pdata*nhaver
269 minhvalid = pdata*nhaver
271 for im in range(self.dataOut.nmodes):
270 for im in range(self.dataOut.nmodes):
272 new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
271 new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
273
272
274
273
275 data_fHeigths_List = []
274 data_fHeigths_List = []
276 data_fZonal_List = []
275 data_fZonal_List = []
277 data_fMeridional_List = []
276 data_fMeridional_List = []
278 data_fVertical_List = []
277 data_fVertical_List = []
279 startDTList = []
278 startDTList = []
280
279
281
280
282 for i in range(ntime):
281 for i in range(ntime):
283 height = old_height
282 height = old_height
284
283
285 start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
284 start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
286 stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
285 stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
287
286
288
287
289 limit_sec1 = time.mktime(start.timetuple())
288 limit_sec1 = time.mktime(start.timetuple())
290 limit_sec2 = time.mktime(stop.timetuple())
289 limit_sec2 = time.mktime(stop.timetuple())
291
290
292 t1 = numpy.where(self.f_timesec >= limit_sec1)
291 t1 = numpy.where(self.f_timesec >= limit_sec1)
293 t2 = numpy.where(self.f_timesec < limit_sec2)
292 t2 = numpy.where(self.f_timesec < limit_sec2)
294 time_select = []
293 time_select = []
295 for val_sec in t1[0]:
294 for val_sec in t1[0]:
296 if val_sec in t2[0]:
295 if val_sec in t2[0]:
297 time_select.append(val_sec)
296 time_select.append(val_sec)
298
297
299
298
300 time_select = numpy.array(time_select,dtype = 'int')
299 time_select = numpy.array(time_select,dtype = 'int')
301 minvalid = numpy.ceil(pdata*nhpoints)
300 minvalid = numpy.ceil(pdata*nhpoints)
302
301
303 zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
302 zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
304 mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
303 mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
305 ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
304 ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
306
305
307 if nhaver > 1:
306 if nhaver > 1:
308 new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
307 new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
309 new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
308 new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
310 new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
309 new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
311
310
312 if len(time_select) > minvalid:
311 if len(time_select) > minvalid:
313 time_average = self.f_timesec[time_select]
312 time_average = self.f_timesec[time_select]
314
313
315 for im in range(self.dataOut.nmodes):
314 for im in range(self.dataOut.nmodes):
316
315
317 for ih in range(self.dataOut.nranges):
316 for ih in range(self.dataOut.nranges):
318 if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
317 if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
319 zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
318 zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
320
319
321 if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
320 if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
322 mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
321 mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
323
322
324 if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
323 if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
325 ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
324 ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
326
325
327 if nhaver > 1:
326 if nhaver > 1:
328 for ih in range(num_hei):
327 for ih in range(num_hei):
329 hvalid = numpy.arange(nhaver) + nhaver*ih
328 hvalid = numpy.arange(nhaver) + nhaver*ih
330
329
331 if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
330 if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
332 new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
331 new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
333
332
334 if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
333 if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
335 new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
334 new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
336
335
337 if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
336 if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
338 new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
337 new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
339 if nhaver > 1:
338 if nhaver > 1:
340 zon_aver = new_zon_aver
339 zon_aver = new_zon_aver
341 mer_aver = new_mer_aver
340 mer_aver = new_mer_aver
342 ver_aver = new_ver_aver
341 ver_aver = new_ver_aver
343 height = new_height
342 height = new_height
344
343
345
344
346 tstart = time_average[0]
345 tstart = time_average[0]
347 tend = time_average[-1]
346 tend = time_average[-1]
348 startTime = time.gmtime(tstart)
347 startTime = time.gmtime(tstart)
349
348
350 year = startTime.tm_year
349 year = startTime.tm_year
351 month = startTime.tm_mon
350 month = startTime.tm_mon
352 day = startTime.tm_mday
351 day = startTime.tm_mday
353 hour = startTime.tm_hour
352 hour = startTime.tm_hour
354 minute = startTime.tm_min
353 minute = startTime.tm_min
355 second = startTime.tm_sec
354 second = startTime.tm_sec
356
355
357 startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
356 startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
358
357
359
358
360 o_height = numpy.array([])
359 o_height = numpy.array([])
361 o_zon_aver = numpy.array([])
360 o_zon_aver = numpy.array([])
362 o_mer_aver = numpy.array([])
361 o_mer_aver = numpy.array([])
363 o_ver_aver = numpy.array([])
362 o_ver_aver = numpy.array([])
364 if self.dataOut.nmodes > 1:
363 if self.dataOut.nmodes > 1:
365 for im in range(self.dataOut.nmodes):
364 for im in range(self.dataOut.nmodes):
366
365
367 if im == 0:
366 if im == 0:
368 h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
367 h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
369 else:
368 else:
370 h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
369 h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
371
370
372
371
373 ht = h_select[0]
372 ht = h_select[0]
374
373
375 o_height = numpy.hstack((o_height,height[im,ht]))
374 o_height = numpy.hstack((o_height,height[im,ht]))
376 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
375 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
377 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
376 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
378 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
377 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
379
378
380 data_fHeigths_List.append(o_height)
379 data_fHeigths_List.append(o_height)
381 data_fZonal_List.append(o_zon_aver)
380 data_fZonal_List.append(o_zon_aver)
382 data_fMeridional_List.append(o_mer_aver)
381 data_fMeridional_List.append(o_mer_aver)
383 data_fVertical_List.append(o_ver_aver)
382 data_fVertical_List.append(o_ver_aver)
384
383
385
384
386 else:
385 else:
387 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
386 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
388 ht = h_select[0]
387 ht = h_select[0]
389 o_height = numpy.hstack((o_height,height[im,ht]))
388 o_height = numpy.hstack((o_height,height[im,ht]))
390 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
389 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
391 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
390 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
392 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
391 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
393
392
394 data_fHeigths_List.append(o_height)
393 data_fHeigths_List.append(o_height)
395 data_fZonal_List.append(o_zon_aver)
394 data_fZonal_List.append(o_zon_aver)
396 data_fMeridional_List.append(o_mer_aver)
395 data_fMeridional_List.append(o_mer_aver)
397 data_fVertical_List.append(o_ver_aver)
396 data_fVertical_List.append(o_ver_aver)
398
397
399
398
400 return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
399 return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
401
400
402
401
General Comments 0
You need to be logged in to leave comments. Login now