##// END OF EJS Templates
jespinoza -
r1361:1dcf28d5b762 merge
parent child
Show More

The requested changes are too big and content was truncated. Show full diff

@@ -1,399 +1,400
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 datetime
9 9 import time
10 10
11 11 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
12 12 from schainpy.model.data.jrodata import Parameters
13 13
14 14
15 15 class BLTRParametersProc(ProcessingUnit):
16 16 '''
17 17 Processing unit for BLTR parameters data (winds)
18 18
19 19 Inputs:
20 20 self.dataOut.nmodes - Number of operation modes
21 21 self.dataOut.nchannels - Number of channels
22 22 self.dataOut.nranges - Number of ranges
23 23
24 24 self.dataOut.data_snr - SNR array
25 25 self.dataOut.data_output - Zonal, Vertical and Meridional velocity array
26 26 self.dataOut.height - Height array (km)
27 27 self.dataOut.time - Time array (seconds)
28 28
29 29 self.dataOut.fileIndex -Index of the file currently read
30 30 self.dataOut.lat - Latitude coordinate of BLTR location
31 31
32 32 self.dataOut.doy - Experiment doy (number of the day in the current year)
33 33 self.dataOut.month - Experiment month
34 34 self.dataOut.day - Experiment day
35 35 self.dataOut.year - Experiment year
36 36 '''
37 37
38 38 def __init__(self):
39 39 '''
40 40 Inputs: None
41 41 '''
42 42 ProcessingUnit.__init__(self)
43 43 self.dataOut = Parameters()
44 44
45 45 def setup(self, mode):
46 46 '''
47 47 '''
48 48 self.dataOut.mode = mode
49 49
50 50 def run(self, mode, snr_threshold=None):
51 51 '''
52 52 Inputs:
53 53 mode = High resolution (0) or Low resolution (1) data
54 54 snr_threshold = snr filter value
55 55 '''
56 56
57 57 if not self.isConfig:
58 58 self.setup(mode)
59 59 self.isConfig = True
60 60
61 61 if self.dataIn.type == 'Parameters':
62 62 self.dataOut.copy(self.dataIn)
63 63
64 64 self.dataOut.data_param = self.dataOut.data[mode]
65 65 self.dataOut.heightList = self.dataOut.height[0]
66 66 self.dataOut.data_snr = self.dataOut.data_snr[mode]
67 67 SNRavg = numpy.average(self.dataOut.data_snr, axis=0)
68 68 SNRavgdB = 10*numpy.log10(SNRavg)
69 69 self.dataOut.data_snr_avg_db = SNRavgdB.reshape(1, *SNRavgdB.shape)
70 70
71 # Censoring Data
71 72 if snr_threshold is not None:
72 73 for i in range(3):
73 74 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
74 75
75 76 # TODO
76 77
77 78 class OutliersFilter(Operation):
78 79
79 80 def __init__(self):
80 81 '''
81 82 '''
82 83 Operation.__init__(self)
83 84
84 85 def run(self, svalue2, method, factor, filter, npoints=9):
85 86 '''
86 87 Inputs:
87 88 svalue - string to select array velocity
88 89 svalue2 - string to choose axis filtering
89 90 method - 0 for SMOOTH or 1 for MEDIAN
90 91 factor - number used to set threshold
91 92 filter - 1 for data filtering using the standard deviation criteria else 0
92 93 npoints - number of points for mask filter
93 94 '''
94 95
95 96 print(' Outliers Filter {} {} / threshold = {}'.format(svalue, svalue, factor))
96 97
97 98
98 99 yaxis = self.dataOut.heightList
99 100 xaxis = numpy.array([[self.dataOut.utctime]])
100 101
101 102 # Zonal
102 103 value_temp = self.dataOut.data_output[0]
103 104
104 105 # Zonal
105 106 value_temp = self.dataOut.data_output[1]
106 107
107 108 # Vertical
108 109 value_temp = numpy.transpose(self.dataOut.data_output[2])
109 110
110 111 htemp = yaxis
111 112 std = value_temp
112 113 for h in range(len(htemp)):
113 114 nvalues_valid = len(numpy.where(numpy.isfinite(value_temp[h]))[0])
114 115 minvalid = npoints
115 116
116 117 #only if valid values greater than the minimum required (10%)
117 118 if nvalues_valid > minvalid:
118 119
119 120 if method == 0:
120 121 #SMOOTH
121 122 w = value_temp[h] - self.Smooth(input=value_temp[h], width=npoints, edge_truncate=1)
122 123
123 124
124 125 if method == 1:
125 126 #MEDIAN
126 127 w = value_temp[h] - self.Median(input=value_temp[h], width = npoints)
127 128
128 129 dw = numpy.std(w[numpy.where(numpy.isfinite(w))],ddof = 1)
129 130
130 131 threshold = dw*factor
131 132 value_temp[numpy.where(w > threshold),h] = numpy.nan
132 133 value_temp[numpy.where(w < -1*threshold),h] = numpy.nan
133 134
134 135
135 136 #At the end
136 137 if svalue2 == 'inHeight':
137 138 value_temp = numpy.transpose(value_temp)
138 139 output_array[:,m] = value_temp
139 140
140 141 if svalue == 'zonal':
141 142 self.dataOut.data_output[0] = output_array
142 143
143 144 elif svalue == 'meridional':
144 145 self.dataOut.data_output[1] = output_array
145 146
146 147 elif svalue == 'vertical':
147 148 self.dataOut.data_output[2] = output_array
148 149
149 150 return self.dataOut.data_output
150 151
151 152
152 153 def Median(self,input,width):
153 154 '''
154 155 Inputs:
155 156 input - Velocity array
156 157 width - Number of points for mask filter
157 158
158 159 '''
159 160
160 161 if numpy.mod(width,2) == 1:
161 162 pc = int((width - 1) / 2)
162 163 cont = 0
163 164 output = []
164 165
165 166 for i in range(len(input)):
166 167 if i >= pc and i < len(input) - pc:
167 168 new2 = input[i-pc:i+pc+1]
168 169 temp = numpy.where(numpy.isfinite(new2))
169 170 new = new2[temp]
170 171 value = numpy.median(new)
171 172 output.append(value)
172 173
173 174 output = numpy.array(output)
174 175 output = numpy.hstack((input[0:pc],output))
175 176 output = numpy.hstack((output,input[-pc:len(input)]))
176 177
177 178 return output
178 179
179 180 def Smooth(self,input,width,edge_truncate = None):
180 181 '''
181 182 Inputs:
182 183 input - Velocity array
183 184 width - Number of points for mask filter
184 185 edge_truncate - 1 for truncate the convolution product else
185 186
186 187 '''
187 188
188 189 if numpy.mod(width,2) == 0:
189 190 real_width = width + 1
190 191 nzeros = width / 2
191 192 else:
192 193 real_width = width
193 194 nzeros = (width - 1) / 2
194 195
195 196 half_width = int(real_width)/2
196 197 length = len(input)
197 198
198 199 gate = numpy.ones(real_width,dtype='float')
199 200 norm_of_gate = numpy.sum(gate)
200 201
201 202 nan_process = 0
202 203 nan_id = numpy.where(numpy.isnan(input))
203 204 if len(nan_id[0]) > 0:
204 205 nan_process = 1
205 206 pb = numpy.zeros(len(input))
206 207 pb[nan_id] = 1.
207 208 input[nan_id] = 0.
208 209
209 210 if edge_truncate == True:
210 211 output = numpy.convolve(input/norm_of_gate,gate,mode='same')
211 212 elif edge_truncate == False or edge_truncate == None:
212 213 output = numpy.convolve(input/norm_of_gate,gate,mode='valid')
213 214 output = numpy.hstack((input[0:half_width],output))
214 215 output = numpy.hstack((output,input[len(input)-half_width:len(input)]))
215 216
216 217 if nan_process:
217 218 pb = numpy.convolve(pb/norm_of_gate,gate,mode='valid')
218 219 pb = numpy.hstack((numpy.zeros(half_width),pb))
219 220 pb = numpy.hstack((pb,numpy.zeros(half_width)))
220 221 output[numpy.where(pb > 0.9999)] = numpy.nan
221 222 input[nan_id] = numpy.nan
222 223 return output
223 224
224 225 def Average(self,aver=0,nhaver=1):
225 226 '''
226 227 Inputs:
227 228 aver - Indicates the time period over which is averaged or consensus data
228 229 nhaver - Indicates the decimation factor in heights
229 230
230 231 '''
231 232 nhpoints = 48
232 233
233 234 lat_piura = -5.17
234 235 lat_huancayo = -12.04
235 236 lat_porcuya = -5.8
236 237
237 238 if '%2.2f'%self.dataOut.lat == '%2.2f'%lat_piura:
238 239 hcm = 3.
239 240 if self.dataOut.year == 2003 :
240 241 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
241 242 nhpoints = 12
242 243
243 244 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_huancayo:
244 245 hcm = 3.
245 246 if self.dataOut.year == 2003 :
246 247 if self.dataOut.doy >= 25 and self.dataOut.doy < 64:
247 248 nhpoints = 12
248 249
249 250
250 251 elif '%2.2f'%self.dataOut.lat == '%2.2f'%lat_porcuya:
251 252 hcm = 5.#2
252 253
253 254 pdata = 0.2
254 255 taver = [1,2,3,4,6,8,12,24]
255 256 t0 = 0
256 257 tf = 24
257 258 ntime =(tf-t0)/taver[aver]
258 259 ti = numpy.arange(ntime)
259 260 tf = numpy.arange(ntime) + taver[aver]
260 261
261 262
262 263 old_height = self.dataOut.heightList
263 264
264 265 if nhaver > 1:
265 266 num_hei = len(self.dataOut.heightList)/nhaver/self.dataOut.nmodes
266 267 deltha = 0.05*nhaver
267 268 minhvalid = pdata*nhaver
268 269 for im in range(self.dataOut.nmodes):
269 270 new_height = numpy.arange(num_hei)*deltha + self.dataOut.height[im,0] + deltha/2.
270 271
271 272
272 273 data_fHeigths_List = []
273 274 data_fZonal_List = []
274 275 data_fMeridional_List = []
275 276 data_fVertical_List = []
276 277 startDTList = []
277 278
278 279
279 280 for i in range(ntime):
280 281 height = old_height
281 282
282 283 start = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(ti[i])) - datetime.timedelta(hours = 5)
283 284 stop = datetime.datetime(self.dataOut.year,self.dataOut.month,self.dataOut.day) + datetime.timedelta(hours = int(tf[i])) - datetime.timedelta(hours = 5)
284 285
285 286
286 287 limit_sec1 = time.mktime(start.timetuple())
287 288 limit_sec2 = time.mktime(stop.timetuple())
288 289
289 290 t1 = numpy.where(self.f_timesec >= limit_sec1)
290 291 t2 = numpy.where(self.f_timesec < limit_sec2)
291 292 time_select = []
292 293 for val_sec in t1[0]:
293 294 if val_sec in t2[0]:
294 295 time_select.append(val_sec)
295 296
296 297
297 298 time_select = numpy.array(time_select,dtype = 'int')
298 299 minvalid = numpy.ceil(pdata*nhpoints)
299 300
300 301 zon_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
301 302 mer_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
302 303 ver_aver = numpy.zeros([self.dataOut.nranges,self.dataOut.nmodes],dtype='f4') + numpy.nan
303 304
304 305 if nhaver > 1:
305 306 new_zon_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
306 307 new_mer_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
307 308 new_ver_aver = numpy.zeros([num_hei,self.dataOut.nmodes],dtype='f4') + numpy.nan
308 309
309 310 if len(time_select) > minvalid:
310 311 time_average = self.f_timesec[time_select]
311 312
312 313 for im in range(self.dataOut.nmodes):
313 314
314 315 for ih in range(self.dataOut.nranges):
315 316 if numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im])) >= minvalid:
316 317 zon_aver[ih,im] = numpy.nansum(self.f_zon[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_zon[time_select,ih,im]))
317 318
318 319 if numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im])) >= minvalid:
319 320 mer_aver[ih,im] = numpy.nansum(self.f_mer[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_mer[time_select,ih,im]))
320 321
321 322 if numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im])) >= minvalid:
322 323 ver_aver[ih,im] = numpy.nansum(self.f_ver[time_select,ih,im]) / numpy.sum(numpy.isfinite(self.f_ver[time_select,ih,im]))
323 324
324 325 if nhaver > 1:
325 326 for ih in range(num_hei):
326 327 hvalid = numpy.arange(nhaver) + nhaver*ih
327 328
328 329 if numpy.sum(numpy.isfinite(zon_aver[hvalid,im])) >= minvalid:
329 330 new_zon_aver[ih,im] = numpy.nansum(zon_aver[hvalid,im]) / numpy.sum(numpy.isfinite(zon_aver[hvalid,im]))
330 331
331 332 if numpy.sum(numpy.isfinite(mer_aver[hvalid,im])) >= minvalid:
332 333 new_mer_aver[ih,im] = numpy.nansum(mer_aver[hvalid,im]) / numpy.sum(numpy.isfinite(mer_aver[hvalid,im]))
333 334
334 335 if numpy.sum(numpy.isfinite(ver_aver[hvalid,im])) >= minvalid:
335 336 new_ver_aver[ih,im] = numpy.nansum(ver_aver[hvalid,im]) / numpy.sum(numpy.isfinite(ver_aver[hvalid,im]))
336 337 if nhaver > 1:
337 338 zon_aver = new_zon_aver
338 339 mer_aver = new_mer_aver
339 340 ver_aver = new_ver_aver
340 341 height = new_height
341 342
342 343
343 344 tstart = time_average[0]
344 345 tend = time_average[-1]
345 346 startTime = time.gmtime(tstart)
346 347
347 348 year = startTime.tm_year
348 349 month = startTime.tm_mon
349 350 day = startTime.tm_mday
350 351 hour = startTime.tm_hour
351 352 minute = startTime.tm_min
352 353 second = startTime.tm_sec
353 354
354 355 startDTList.append(datetime.datetime(year,month,day,hour,minute,second))
355 356
356 357
357 358 o_height = numpy.array([])
358 359 o_zon_aver = numpy.array([])
359 360 o_mer_aver = numpy.array([])
360 361 o_ver_aver = numpy.array([])
361 362 if self.dataOut.nmodes > 1:
362 363 for im in range(self.dataOut.nmodes):
363 364
364 365 if im == 0:
365 366 h_select = numpy.where(numpy.bitwise_and(height[0,:] >=0,height[0,:] <= hcm,numpy.isfinite(height[0,:])))
366 367 else:
367 368 h_select = numpy.where(numpy.bitwise_and(height[1,:] > hcm,height[1,:] < 20,numpy.isfinite(height[1,:])))
368 369
369 370
370 371 ht = h_select[0]
371 372
372 373 o_height = numpy.hstack((o_height,height[im,ht]))
373 374 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
374 375 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
375 376 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
376 377
377 378 data_fHeigths_List.append(o_height)
378 379 data_fZonal_List.append(o_zon_aver)
379 380 data_fMeridional_List.append(o_mer_aver)
380 381 data_fVertical_List.append(o_ver_aver)
381 382
382 383
383 384 else:
384 385 h_select = numpy.where(numpy.bitwise_and(height[0,:] <= hcm,numpy.isfinite(height[0,:])))
385 386 ht = h_select[0]
386 387 o_height = numpy.hstack((o_height,height[im,ht]))
387 388 o_zon_aver = numpy.hstack((o_zon_aver,zon_aver[ht,im]))
388 389 o_mer_aver = numpy.hstack((o_mer_aver,mer_aver[ht,im]))
389 390 o_ver_aver = numpy.hstack((o_ver_aver,ver_aver[ht,im]))
390 391
391 392 data_fHeigths_List.append(o_height)
392 393 data_fZonal_List.append(o_zon_aver)
393 394 data_fMeridional_List.append(o_mer_aver)
394 395 data_fVertical_List.append(o_ver_aver)
395 396
396 397
397 398 return startDTList, data_fHeigths_List, data_fZonal_List, data_fMeridional_List, data_fVertical_List
398 399
399 400
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now