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