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