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