##// END OF EJS Templates
schain
RSS
Merge branch 'master' of http://jro-dev.igp.gob.pe/rhodecode/schain
jespinoza -
r1361:1dcf28d5b762 merge
parent child
Show More
Side by Side Unified
Context file:
r1361:1dcf28d5b762 Loading diff...:
Show file before | Show file after | Hide comments
@@ -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
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1164 lines changed) Show them Hide them
@@ -1,4000 +1,3886
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import scipy
5 5 import re
6 6 import datetime
7 7 import copy
8 8 import sys
9 9 import importlib
10 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 13 import time
14 14
15 15 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
16 16 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
17 17 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
18 18 from scipy import asarray as ar,exp
19 19 from scipy.optimize import curve_fit
20 20 from schainpy.utils import log
21 21 import warnings
22 22 from numpy import NaN
23 23 from scipy.optimize.optimize import OptimizeWarning
24 24 warnings.filterwarnings('ignore')
25 25
26 26 import matplotlib.pyplot as plt
27 27
28 28 SPEED_OF_LIGHT = 299792458
29 29
30
31 30 '''solving pickling issue'''
32 31
33 32 def _pickle_method(method):
34 33 func_name = method.__func__.__name__
35 34 obj = method.__self__
36 35 cls = method.__self__.__class__
37 36 return _unpickle_method, (func_name, obj, cls)
38 37
39 38 def _unpickle_method(func_name, obj, cls):
40 39 for cls in cls.mro():
41 40 try:
42 41 func = cls.__dict__[func_name]
43 42 except KeyError:
44 43 pass
45 44 else:
46 45 break
47 46 return func.__get__(obj, cls)
48 47
49 48
50 49 class ParametersProc(ProcessingUnit):
51 50
52 51 METHODS = {}
53 52 nSeconds = None
54 53
55 54 def __init__(self):
56 55 ProcessingUnit.__init__(self)
57 56
58 57 # self.objectDict = {}
59 58 self.buffer = None
60 59 self.firstdatatime = None
61 60 self.profIndex = 0
62 61 self.dataOut = Parameters()
63 62 self.setupReq = False #Agregar a todas las unidades de proc
64 63
65 64 def __updateObjFromInput(self):
66 65
67 66 self.dataOut.inputUnit = self.dataIn.type
68 67
69 68 self.dataOut.timeZone = self.dataIn.timeZone
70 69 self.dataOut.dstFlag = self.dataIn.dstFlag
71 70 self.dataOut.errorCount = self.dataIn.errorCount
72 71 self.dataOut.useLocalTime = self.dataIn.useLocalTime
73 72
74 73 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
75 74 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
76 75 self.dataOut.channelList = self.dataIn.channelList
77 76 self.dataOut.heightList = self.dataIn.heightList
78 77 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
79 78 # self.dataOut.nHeights = self.dataIn.nHeights
80 79 # self.dataOut.nChannels = self.dataIn.nChannels
81 80 # self.dataOut.nBaud = self.dataIn.nBaud
82 81 # self.dataOut.nCode = self.dataIn.nCode
83 82 # self.dataOut.code = self.dataIn.code
84 83 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
85 84 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
86 85 # self.dataOut.utctime = self.firstdatatime
87 86 self.dataOut.utctime = self.dataIn.utctime
88 87 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
89 88 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
90 89 self.dataOut.nCohInt = self.dataIn.nCohInt
91 90 # self.dataOut.nIncohInt = 1
92 91 # self.dataOut.ippSeconds = self.dataIn.ippSeconds
93 92 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
94 93 self.dataOut.timeInterval1 = self.dataIn.timeInterval
95 94 self.dataOut.heightList = self.dataIn.heightList
96 95 self.dataOut.frequency = self.dataIn.frequency
97 96 # self.dataOut.noise = self.dataIn.noise
98 97
99 98 def run(self):
100 99
101 100
102 101
103 102 #---------------------- Voltage Data ---------------------------
104 103
105 104 if self.dataIn.type == "Voltage":
106 105
107 106 self.__updateObjFromInput()
108 107 self.dataOut.data_pre = self.dataIn.data.copy()
109 108 self.dataOut.flagNoData = False
110 109 self.dataOut.utctimeInit = self.dataIn.utctime
111 110 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
112 111 if hasattr(self.dataIn, 'dataPP_POW'):
113 112 self.dataOut.dataPP_POW = self.dataIn.dataPP_POW
114 113
115 114 if hasattr(self.dataIn, 'dataPP_POWER'):
116 115 self.dataOut.dataPP_POWER = self.dataIn.dataPP_POWER
117 116
118 117 if hasattr(self.dataIn, 'dataPP_DOP'):
119 118 self.dataOut.dataPP_DOP = self.dataIn.dataPP_DOP
120 119
121 120 if hasattr(self.dataIn, 'dataPP_SNR'):
122 121 self.dataOut.dataPP_SNR = self.dataIn.dataPP_SNR
123 122
124 123 if hasattr(self.dataIn, 'dataPP_WIDTH'):
125 124 self.dataOut.dataPP_WIDTH = self.dataIn.dataPP_WIDTH
126 125 return
127 126
128 127 #---------------------- Spectra Data ---------------------------
129 128
130 129 if self.dataIn.type == "Spectra":
131 130
132 self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
131 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
133 132 self.dataOut.data_spc = self.dataIn.data_spc
134 133 self.dataOut.data_cspc = self.dataIn.data_cspc
135 134 self.dataOut.nProfiles = self.dataIn.nProfiles
136 135 self.dataOut.nIncohInt = self.dataIn.nIncohInt
137 136 self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
138 137 self.dataOut.ippFactor = self.dataIn.ippFactor
139 138 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
140 139 self.dataOut.spc_noise = self.dataIn.getNoise()
141 140 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
142 141 # self.dataOut.normFactor = self.dataIn.normFactor
143 142 self.dataOut.pairsList = self.dataIn.pairsList
144 143 self.dataOut.groupList = self.dataIn.pairsList
145 144 self.dataOut.flagNoData = False
146 145
147 146 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
148 147 self.dataOut.ChanDist = self.dataIn.ChanDist
149 148 else: self.dataOut.ChanDist = None
150 149
151 150 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
152 151 # self.dataOut.VelRange = self.dataIn.VelRange
153 152 #else: self.dataOut.VelRange = None
154 153
155 154 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
156 155 self.dataOut.RadarConst = self.dataIn.RadarConst
157 156
158 157 if hasattr(self.dataIn, 'NPW'): #NPW
159 158 self.dataOut.NPW = self.dataIn.NPW
160 159
161 160 if hasattr(self.dataIn, 'COFA'): #COFA
162 161 self.dataOut.COFA = self.dataIn.COFA
163 162
164 163
165 164
166 165 #---------------------- Correlation Data ---------------------------
167 166
168 167 if self.dataIn.type == "Correlation":
169 168 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
170 169
171 170 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
172 171 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
173 172 self.dataOut.groupList = (acf_pairs, ccf_pairs)
174 173
175 174 self.dataOut.abscissaList = self.dataIn.lagRange
176 175 self.dataOut.noise = self.dataIn.noise
177 176 self.dataOut.data_snr = self.dataIn.SNR
178 177 self.dataOut.flagNoData = False
179 178 self.dataOut.nAvg = self.dataIn.nAvg
180 179
181 180 #---------------------- Parameters Data ---------------------------
182 181
183 182 if self.dataIn.type == "Parameters":
184 183 self.dataOut.copy(self.dataIn)
185 184 self.dataOut.flagNoData = False
186 185
187 186 return True
188 187
189 188 self.__updateObjFromInput()
190 189 self.dataOut.utctimeInit = self.dataIn.utctime
191 190 self.dataOut.paramInterval = self.dataIn.timeInterval
192 191
193 192 return
194 193
195 194
196 195 def target(tups):
197 196
198 197 obj, args = tups
199 198
200 199 return obj.FitGau(args)
201 200
201 class RemoveWideGC(Operation):
202 ''' This class remove the wide clutter and replace it with a simple interpolation points
203 This mainly applies to CLAIRE radar
202 204
203 class SpectralFilters(Operation):
204
205 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
206
207 LimitR : It is the limit in m/s of Rainfall
208 LimitW : It is the limit in m/s for Winds
205 ClutterWidth : Width to look for the clutter peak
209 206
210 207 Input:
211 208
212 209 self.dataOut.data_pre : SPC and CSPC
213 210 self.dataOut.spc_range : To select wind and rainfall velocities
214 211
215 212 Affected:
216 213
217 214 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
218 self.dataOut.spcparam_range : Used in SpcParamPlot
219 self.dataOut.SPCparam : Used in PrecipitationProc
220
221 215
216 Written by D. ScipiΓ³n 25.02.2021
222 217 '''
223
224 218 def __init__(self):
225 219 Operation.__init__(self)
226 self.i=0
227
228 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
229
230
231 #Limite de vientos
232 LimitR = PositiveLimit
233 LimitN = NegativeLimit
220 self.i = 0
221 self.ich = 0
222 self.ir = 0
223
224 def run(self, dataOut, ClutterWidth=2.5):
225 # print ('Entering RemoveWideGC ... ')
234 226
235 227 self.spc = dataOut.data_pre[0].copy()
236 self.cspc = dataOut.data_pre[1].copy()
237
238 self.Num_Hei = self.spc.shape[2]
239 self.Num_Bin = self.spc.shape[1]
228 self.spc_out = dataOut.data_pre[0].copy()
240 229 self.Num_Chn = self.spc.shape[0]
230 self.Num_Hei = self.spc.shape[2]
231 VelRange = dataOut.spc_range[2][:-1]
232 dv = VelRange[1]-VelRange[0]
233
234 # Find the velocities that corresponds to zero
235 gc_values = numpy.squeeze(numpy.where(numpy.abs(VelRange) <= ClutterWidth))
236
237 # Removing novalid data from the spectra
238 for ich in range(self.Num_Chn) :
239 for ir in range(self.Num_Hei) :
240 # Estimate the noise at each range
241 HSn = hildebrand_sekhon(self.spc[ich,:,ir],dataOut.nIncohInt)
242
243 # Removing the noise floor at each range
244 novalid = numpy.where(self.spc[ich,:,ir] < HSn)
245 self.spc[ich,novalid,ir] = HSn
246
247 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
248 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
249 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
250 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
251 continue
252 junk3 = numpy.squeeze(numpy.diff(j1index))
253 junk4 = numpy.squeeze(numpy.diff(j2index))
254
255 valleyindex = j2index[numpy.where(junk4>1)]
256 peakindex = j1index[numpy.where(junk3>1)]
257
258 isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
259 if numpy.size(isvalid) == 0 :
260 continue
261 if numpy.size(isvalid) >1 :
262 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
263 isvalid = isvalid[vindex]
264
265 # clutter peak
266 gcpeak = peakindex[isvalid]
267 vl = numpy.where(valleyindex < gcpeak)
268 if numpy.size(vl) == 0:
269 continue
270 gcvl = valleyindex[vl[0][-1]]
271 vr = numpy.where(valleyindex > gcpeak)
272 if numpy.size(vr) == 0:
273 continue
274 gcvr = valleyindex[vr[0][0]]
275
276 # Removing the clutter
277 interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
278 gcindex = gc_values[gcvl+1:gcvr-1]
279 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
280
281 dataOut.data_pre[0] = self.spc_out
282 #print ('Leaving RemoveWideGC ... ')
283 return dataOut
241 284
242 VelRange = dataOut.spc_range[2]
243 TimeRange = dataOut.spc_range[1]
244 FrecRange = dataOut.spc_range[0]
245
246 Vmax= 2*numpy.max(dataOut.spc_range[2])
247 Tmax= 2*numpy.max(dataOut.spc_range[1])
248 Fmax= 2*numpy.max(dataOut.spc_range[0])
249
250 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
251 Breaker1R=numpy.where(VelRange == Breaker1R)
252
253 Delta = self.Num_Bin/2 - Breaker1R[0]
254
255
256 '''Reacomodando SPCrange'''
257
258 VelRange=numpy.roll(VelRange,-(int(self.Num_Bin/2)) ,axis=0)
259
260 VelRange[-(int(self.Num_Bin/2)):]+= Vmax
285 class SpectralFilters(Operation):
286 ''' This class allows to replace the novalid values with noise for each channel
287 This applies to CLAIRE RADAR
261 288
262 FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0)
289 PositiveLimit : RightLimit of novalid data
290 NegativeLimit : LeftLimit of novalid data
263 291
264 FrecRange[-(int(self.Num_Bin/2)):]+= Fmax
292 Input:
265 293
266 TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0)
294 self.dataOut.data_pre : SPC and CSPC
295 self.dataOut.spc_range : To select wind and rainfall velocities
267 296
268 TimeRange[-(int(self.Num_Bin/2)):]+= Tmax
297 Affected:
269 298
270 ''' ------------------ '''
299 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
271 300
272 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
273 Breaker2R=numpy.where(VelRange == Breaker2R)
301 Written by D. ScipiΓ³n 29.01.2021
302 '''
303 def __init__(self):
304 Operation.__init__(self)
305 self.i = 0
306
307 def run(self, dataOut, ):
274 308
309 self.spc = dataOut.data_pre[0].copy()
310 self.Num_Chn = self.spc.shape[0]
311 VelRange = dataOut.spc_range[2]
275 312
276 SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1)
313 # novalid corresponds to data within the Negative and PositiveLimit
314
277 315
278 SPCcut = SPCroll.copy()
316 # Removing novalid data from the spectra
279 317 for i in range(self.Num_Chn):
280
281 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
282 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
283
284 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
285 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
286
287 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
288 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
289
290 SPC_ch1 = SPCroll
291
292 SPC_ch2 = SPCcut
293
294 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
295 dataOut.SPCparam = numpy.asarray(SPCparam)
296
297
298 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
299
300 dataOut.spcparam_range[2]=VelRange
301 dataOut.spcparam_range[1]=TimeRange
302 dataOut.spcparam_range[0]=FrecRange
318 self.spc[i,novalid,:] = dataOut.noise[i]
319 dataOut.data_pre[0] = self.spc
303 320 return dataOut
304 321
305 322 class GaussianFit(Operation):
306 323
307 324 '''
308 325 Function that fit of one and two generalized gaussians (gg) based
309 326 on the PSD shape across an "power band" identified from a cumsum of
310 327 the measured spectrum - noise.
311 328
312 329 Input:
313 330 self.dataOut.data_pre : SelfSpectra
314 331
315 332 Output:
316 333 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
317 334
318 335 '''
319 336 def __init__(self):
320 337 Operation.__init__(self)
321 338 self.i=0
322 339
323 340
324 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
341 # def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
342 def run(self, dataOut, SNRdBlimit=-9, method='generalized'):
325 343 """This routine will find a couple of generalized Gaussians to a power spectrum
344 methods: generalized, squared
326 345 input: spc
327 346 output:
328 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
347 noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
329 348 """
330
349 print ('Entering ',method,' double Gaussian fit')
331 350 self.spc = dataOut.data_pre[0].copy()
332 351 self.Num_Hei = self.spc.shape[2]
333 352 self.Num_Bin = self.spc.shape[1]
334 353 self.Num_Chn = self.spc.shape[0]
335 Vrange = dataOut.abscissaList
336
337 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
338 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
339 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
340 SPC_ch1[:] = numpy.NaN
341 SPC_ch2[:] = numpy.NaN
342
343 354
344 355 start_time = time.time()
345 356
346 noise_ = dataOut.spc_noise[0].copy()
347
348
349 357 pool = Pool(processes=self.Num_Chn)
350 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
358 args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
351 359 objs = [self for __ in range(self.Num_Chn)]
352 360 attrs = list(zip(objs, args))
353 gauSPC = pool.map(target, attrs)
354 dataOut.SPCparam = numpy.asarray(SPCparam)
355
356 ''' Parameters:
357 1. Amplitude
358 2. Shift
359 3. Width
360 4. Power
361 '''
361 DGauFitParam = pool.map(target, attrs)
362 # Parameters:
363 # 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
364 dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
365
366 # Double Gaussian Curves
367 gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
368 gau0[:] = numpy.NaN
369 gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
370 gau1[:] = numpy.NaN
371 x_mtr = numpy.transpose(numpy.tile(dataOut.getVelRange(1)[:-1], (self.Num_Hei,1)))
372 for iCh in range(self.Num_Chn):
373 N0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,0]] * self.Num_Bin))
374 N1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,1]] * self.Num_Bin))
375 A0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,0]] * self.Num_Bin))
376 A1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,1]] * self.Num_Bin))
377 v0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,0]] * self.Num_Bin))
378 v1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,1]] * self.Num_Bin))
379 s0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,0]] * self.Num_Bin))
380 s1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,1]] * self.Num_Bin))
381 if method == 'genealized':
382 p0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,0]] * self.Num_Bin))
383 p1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,1]] * self.Num_Bin))
384 elif method == 'squared':
385 p0 = 2.
386 p1 = 2.
387 gau0[iCh] = A0*numpy.exp(-0.5*numpy.abs((x_mtr-v0)/s0)**p0)+N0
388 gau1[iCh] = A1*numpy.exp(-0.5*numpy.abs((x_mtr-v1)/s1)**p1)+N1
389 dataOut.GaussFit0 = gau0
390 dataOut.GaussFit1 = gau1
391
392 print('Leaving ',method ,' double Gaussian fit')
393 return dataOut
362 394
363 395 def FitGau(self, X):
364
365 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
366
367 SPCparam = []
368 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
369 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
370 SPC_ch1[:] = 0#numpy.NaN
371 SPC_ch2[:] = 0#numpy.NaN
372
373
374
396 # print('Entering FitGau')
397 # Assigning the variables
398 Vrange, ch, wnoise, num_intg, SNRlimit = X
399 # Noise Limits
400 noisebl = wnoise * 0.9
401 noisebh = wnoise * 1.1
402 # Radar Velocity
403 Va = max(Vrange)
404 deltav = Vrange[1] - Vrange[0]
405 x = numpy.arange(self.Num_Bin)
406
407 # print ('stop 0')
408
409 # 5 parameters, 2 Gaussians
410 DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
411 DGauFitParam[:] = numpy.NaN
412
413 # SPCparam = []
414 # SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
415 # SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
416 # SPC_ch1[:] = 0 #numpy.NaN
417 # SPC_ch2[:] = 0 #numpy.NaN
418 # print ('stop 1')
375 419 for ht in range(self.Num_Hei):
376
377
420 # print (ht)
421 # print ('stop 2')
422 # Spectra at each range
378 423 spc = numpy.asarray(self.spc)[ch,:,ht]
424 snr = ( spc.mean() - wnoise ) / wnoise
425 snrdB = 10.*numpy.log10(snr)
379 426
427 #print ('stop 3')
428 if snrdB < SNRlimit :
429 # snr = numpy.NaN
430 # SPC_ch1[:,ht] = 0#numpy.NaN
431 # SPC_ch1[:,ht] = 0#numpy.NaN
432 # SPCparam = (SPC_ch1,SPC_ch2)
433 # print ('SNR less than SNRth')
434 continue
435 # wnoise = hildebrand_sekhon(spc,num_intg)
436 # print ('stop 2.01')
380 437 #############################################
381 438 # normalizing spc and noise
382 439 # This part differs from gg1
383 spc_norm_max = max(spc)
440 # spc_norm_max = max(spc) #commented by D. ScipiΓ³n 19.03.2021
384 441 #spc = spc / spc_norm_max
385 pnoise = pnoise #/ spc_norm_max
442 # pnoise = pnoise #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
386 443 #############################################
387 444
445 # print ('stop 2.1')
388 446 fatspectra=1.0
389
390 wnoise = noise_ #/ spc_norm_max
447 # noise per channel.... we might want to use the noise at each range
448
449 # wnoise = noise_ #/ spc_norm_max #commented by D. ScipiΓ³n 19.03.2021
391 450 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
392 451 #if wnoise>1.1*pnoise: # to be tested later
393 452 # wnoise=pnoise
394 noisebl=wnoise*0.9;
395 noisebh=wnoise*1.1
396 spc=spc-wnoise
453 # noisebl = wnoise*0.9
454 # noisebh = wnoise*1.1
455 spc = spc - wnoise # signal
397 456
398 minx=numpy.argmin(spc)
457 # print ('stop 2.2')
458 minx = numpy.argmin(spc)
399 459 #spcs=spc.copy()
400 spcs=numpy.roll(spc,-minx)
401 cum=numpy.cumsum(spcs)
402 tot_noise=wnoise * self.Num_Bin #64;
403
404 snr = sum(spcs)/tot_noise
405 snrdB=10.*numpy.log10(snr)
406
407 if snrdB < SNRlimit :
408 snr = numpy.NaN
409 SPC_ch1[:,ht] = 0#numpy.NaN
410 SPC_ch1[:,ht] = 0#numpy.NaN
411 SPCparam = (SPC_ch1,SPC_ch2)
412 continue
460 spcs = numpy.roll(spc,-minx)
461 cum = numpy.cumsum(spcs)
462 # tot_noise = wnoise * self.Num_Bin #64;
463
464 # print ('stop 2.3')
465 # snr = sum(spcs) / tot_noise
466 # snrdB = 10.*numpy.log10(snr)
467 #print ('stop 3')
468 # if snrdB < SNRlimit :
469 # snr = numpy.NaN
470 # SPC_ch1[:,ht] = 0#numpy.NaN
471 # SPC_ch1[:,ht] = 0#numpy.NaN
472 # SPCparam = (SPC_ch1,SPC_ch2)
473 # print ('SNR less than SNRth')
474 # continue
413 475
414 476
415 477 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
416 478 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
417
418 cummax=max(cum);
419 epsi=0.08*fatspectra # cumsum to narrow down the energy region
420 cumlo=cummax*epsi;
421 cumhi=cummax*(1-epsi)
422 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
423
424
479 # print ('stop 4')
480 cummax = max(cum)
481 epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
482 cumlo = cummax * epsi
483 cumhi = cummax * (1-epsi)
484 powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
485
486 # print ('stop 5')
425 487 if len(powerindex) < 1:# case for powerindex 0
488 # print ('powerindex < 1')
426 489 continue
427 powerlo=powerindex[0]
428 powerhi=powerindex[-1]
429 powerwidth=powerhi-powerlo
430
431 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
432 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
433 midpeak=(firstpeak+secondpeak)/2.
434 firstamp=spcs[int(firstpeak)]
435 secondamp=spcs[int(secondpeak)]
436 midamp=spcs[int(midpeak)]
490 powerlo = powerindex[0]
491 powerhi = powerindex[-1]
492 powerwidth = powerhi-powerlo
493 if powerwidth <= 1:
494 # print('powerwidth <= 1')
495 continue
496
497 # print ('stop 6')
498 firstpeak = powerlo + powerwidth/10.# first gaussian energy location
499 secondpeak = powerhi - powerwidth/10. #second gaussian energy location
500 midpeak = (firstpeak + secondpeak)/2.
501 firstamp = spcs[int(firstpeak)]
502 secondamp = spcs[int(secondpeak)]
503 midamp = spcs[int(midpeak)]
437 504
438 x=numpy.arange( self.Num_Bin )
439 y_data=spc+wnoise
505 y_data = spc + wnoise
440 506
441 507 ''' single Gaussian '''
442 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
443 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
444 power0=2.
445 amplitude0=midamp
446 state0=[shift0,width0,amplitude0,power0,wnoise]
447 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
448 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
449
450 chiSq1=lsq1[1];
451
452
508 shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
509 width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
510 power0 = 2.
511 amplitude0 = midamp
512 state0 = [shift0,width0,amplitude0,power0,wnoise]
513 bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
514 lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
515 # print ('stop 7.1')
516 # print (bnds)
517
518 chiSq1=lsq1[1]
519
520 # print ('stop 8')
453 521 if fatspectra<1.0 and powerwidth<4:
454 522 choice=0
455 523 Amplitude0=lsq1[0][2]
456 524 shift0=lsq1[0][0]
457 525 width0=lsq1[0][1]
458 526 p0=lsq1[0][3]
459 527 Amplitude1=0.
460 528 shift1=0.
461 529 width1=0.
462 530 p1=0.
463 531 noise=lsq1[0][4]
464 532 #return (numpy.array([shift0,width0,Amplitude0,p0]),
465 533 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
466
467 ''' two gaussians '''
534
535 # print ('stop 9')
536 ''' two Gaussians '''
468 537 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
469 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
470 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
471 width0=powerwidth/6.;
472 width1=width0
473 power0=2.;
474 power1=power0
475 amplitude0=firstamp;
476 amplitude1=secondamp
477 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
538 shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
539 shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
540 width0 = powerwidth/6.
541 width1 = width0
542 power0 = 2.
543 power1 = power0
544 amplitude0 = firstamp
545 amplitude1 = secondamp
546 state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
478 547 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
479 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
548 bnds=((0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
480 549 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
481 550
551 # print ('stop 10')
482 552 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
483 553
554 # print ('stop 11')
555 chiSq2 = lsq2[1]
484 556
485 chiSq2=lsq2[1];
486
487
557 # print ('stop 12')
488 558
489 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
559 oneG = (chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
490 560
561 # print ('stop 13')
491 562 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
492 563 if oneG:
493 choice=0
564 choice = 0
494 565 else:
495 w1=lsq2[0][1]; w2=lsq2[0][5]
496 a1=lsq2[0][2]; a2=lsq2[0][6]
497 p1=lsq2[0][3]; p2=lsq2[0][7]
498 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
499 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
500 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
566 w1 = lsq2[0][1]; w2 = lsq2[0][5]
567 a1 = lsq2[0][2]; a2 = lsq2[0][6]
568 p1 = lsq2[0][3]; p2 = lsq2[0][7]
569 s1 = (2**(1+1./p1))*scipy.special.gamma(1./p1)/p1
570 s2 = (2**(1+1./p2))*scipy.special.gamma(1./p2)/p2
571 gp1 = a1*w1*s1; gp2 = a2*w2*s2 # power content of each ggaussian with proper p scaling
501 572
502 573 if gp1>gp2:
503 574 if a1>0.7*a2:
504 choice=1
575 choice = 1
505 576 else:
506 choice=2
577 choice = 2
507 578 elif gp2>gp1:
508 579 if a2>0.7*a1:
509 choice=2
580 choice = 2
510 581 else:
511 choice=1
582 choice = 1
512 583 else:
513 choice=numpy.argmax([a1,a2])+1
584 choice = numpy.argmax([a1,a2])+1
514 585 #else:
515 586 #choice=argmin([std2a,std2b])+1
516 587
517 588 else: # with low SNR go to the most energetic peak
518 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
519
520
521 shift0=lsq2[0][0];
522 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
523 shift1=lsq2[0][4];
524 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
525
526 max_vel = 1.0
527
589 choice = numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
590
591 # print ('stop 14')
592 shift0 = lsq2[0][0]
593 vel0 = Vrange[0] + shift0 * deltav
594 shift1 = lsq2[0][4]
595 # vel1=Vrange[0] + shift1 * deltav
596
597 # max_vel = 1.0
598 # Va = max(Vrange)
599 # deltav = Vrange[1]-Vrange[0]
600 # print ('stop 15')
528 601 #first peak will be 0, second peak will be 1
529 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
530 shift0=lsq2[0][0]
531 width0=lsq2[0][1]
532 Amplitude0=lsq2[0][2]
533 p0=lsq2[0][3]
534
535 shift1=lsq2[0][4]
536 width1=lsq2[0][5]
537 Amplitude1=lsq2[0][6]
538 p1=lsq2[0][7]
539 noise=lsq2[0][8]
602 # if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.ScipiΓ³n 19.03.2021
603 if vel0 > -Va and vel0 < Va : #first peak is in the correct range
604 shift0 = lsq2[0][0]
605 width0 = lsq2[0][1]
606 Amplitude0 = lsq2[0][2]
607 p0 = lsq2[0][3]
608
609 shift1 = lsq2[0][4]
610 width1 = lsq2[0][5]
611 Amplitude1 = lsq2[0][6]
612 p1 = lsq2[0][7]
613 noise = lsq2[0][8]
540 614 else:
541 shift1=lsq2[0][0]
542 width1=lsq2[0][1]
543 Amplitude1=lsq2[0][2]
544 p1=lsq2[0][3]
615 shift1 = lsq2[0][0]
616 width1 = lsq2[0][1]
617 Amplitude1 = lsq2[0][2]
618 p1 = lsq2[0][3]
545 619
546 shift0=lsq2[0][4]
547 width0=lsq2[0][5]
548 Amplitude0=lsq2[0][6]
549 p0=lsq2[0][7]
550 noise=lsq2[0][8]
620 shift0 = lsq2[0][4]
621 width0 = lsq2[0][5]
622 Amplitude0 = lsq2[0][6]
623 p0 = lsq2[0][7]
624 noise = lsq2[0][8]
551 625
552 626 if Amplitude0<0.05: # in case the peak is noise
553 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
627 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
554 628 if Amplitude1<0.05:
555 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
556
557
558 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
559 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
560 SPCparam = (SPC_ch1,SPC_ch2)
561
562
563 return GauSPC
629 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
630
631 # print ('stop 16 ')
632 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
633 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
634 # SPCparam = (SPC_ch1,SPC_ch2)
635
636 DGauFitParam[0,ht,0] = noise
637 DGauFitParam[0,ht,1] = noise
638 DGauFitParam[1,ht,0] = Amplitude0
639 DGauFitParam[1,ht,1] = Amplitude1
640 DGauFitParam[2,ht,0] = Vrange[0] + shift0 * deltav
641 DGauFitParam[2,ht,1] = Vrange[0] + shift1 * deltav
642 DGauFitParam[3,ht,0] = width0 * deltav
643 DGauFitParam[3,ht,1] = width1 * deltav
644 DGauFitParam[4,ht,0] = p0
645 DGauFitParam[4,ht,1] = p1
646
647 # print (DGauFitParam.shape)
648 # print ('Leaving FitGau')
649 return DGauFitParam
650 # return SPCparam
651 # return GauSPC
564 652
565 653 def y_model1(self,x,state):
566 shift0,width0,amplitude0,power0,noise=state
567 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
568
569 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
570
571 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
572 return model0+model0u+model0d+noise
654 shift0, width0, amplitude0, power0, noise = state
655 model0 = amplitude0*numpy.exp(-0.5*abs((x - shift0)/width0)**power0)
656 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
657 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
658 return model0 + model0u + model0d + noise
573 659
574 660 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
575 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
576 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
577
578 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
579
580 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
581 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
582
583 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
584
585 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
586 return model0+model0u+model0d+model1+model1u+model1d+noise
661 shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
662 model0 = amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
663 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
664 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
665
666 model1 = amplitude1*numpy.exp(-0.5*abs((x - shift1)/width1)**power1)
667 model1u = amplitude1*numpy.exp(-0.5*abs((x - shift1 - self.Num_Bin)/width1)**power1)
668 model1d = amplitude1*numpy.exp(-0.5*abs((x - shift1 + self.Num_Bin)/width1)**power1)
669 return model0 + model0u + model0d + model1 + model1u + model1d + noise
587 670
588 671 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
589 672
590 673 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
591 674
592 675 def misfit2(self,state,y_data,x,num_intg):
593 676 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
594 677
595 678
596 679
597 680 class PrecipitationProc(Operation):
598 681
599 682 '''
600 683 Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
601 684
602 685 Input:
603 686 self.dataOut.data_pre : SelfSpectra
604 687
605 688 Output:
606 689
607 690 self.dataOut.data_output : Reflectivity factor, rainfall Rate
608 691
609 692
610 693 Parameters affected:
611 694 '''
612 695
613 696 def __init__(self):
614 697 Operation.__init__(self)
615 698 self.i=0
616 699
617
618 def gaus(self,xSamples,Amp,Mu,Sigma):
619 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
620
621
622
623 def Moments(self, ySamples, xSamples):
624 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
625 yNorm = ySamples / Pot
626
627 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
628 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
629 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
630
631 return numpy.array([Pot, Vr, Desv])
632
633 700 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
634 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
635
701 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30):
636 702
637 Velrange = dataOut.spcparam_range[2]
638 FrecRange = dataOut.spcparam_range[0]
639
640 dV= Velrange[1]-Velrange[0]
641 dF= FrecRange[1]-FrecRange[0]
703 # print ('Entering PrecepitationProc ... ')
642 704
643 705 if radar == "MIRA35C" :
644 706
645 707 self.spc = dataOut.data_pre[0].copy()
646 708 self.Num_Hei = self.spc.shape[2]
647 709 self.Num_Bin = self.spc.shape[1]
648 710 self.Num_Chn = self.spc.shape[0]
649 711 Ze = self.dBZeMODE2(dataOut)
650 712
651 713 else:
652 714
653 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
654
655 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
715 self.spc = dataOut.data_pre[0].copy()
656 716
717 #NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
657 718 self.spc[:,:,0:7]= numpy.NaN
658 719
659 """##########################################"""
660
661 720 self.Num_Hei = self.spc.shape[2]
662 721 self.Num_Bin = self.spc.shape[1]
663 722 self.Num_Chn = self.spc.shape[0]
664 723
724 VelRange = dataOut.spc_range[2]
725
665 726 ''' Se obtiene la constante del RADAR '''
666 727
667 728 self.Pt = Pt
668 729 self.Gt = Gt
669 730 self.Gr = Gr
670 731 self.Lambda = Lambda
671 732 self.aL = aL
672 733 self.tauW = tauW
673 self.ThetaT = ThetaT
734 self.ThetaT = ThetaT
674 735 self.ThetaR = ThetaR
736 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
737 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
738 self.lr = 10**(5.73/10) # Perdida en cables Rx 5.73 dB
675 739
676 740 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
677 741 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
678 742 RadarConstant = 10e-26 * Numerator / Denominator #
679
680 ''' ============================= '''
681
682 self.spc[0] = (self.spc[0]-dataOut.noise[0])
683 self.spc[1] = (self.spc[1]-dataOut.noise[1])
684 self.spc[2] = (self.spc[2]-dataOut.noise[2])
685
686 self.spc[ numpy.where(self.spc < 0)] = 0
687
688 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
689 SPCmean[ numpy.where(SPCmean < 0)] = 0
690
691 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
692 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
693 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
694
695 Pr = SPCmean[:,:]
696
697 VelMeteoro = numpy.mean(SPCmean,axis=0)
698
699 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
700 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
701 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
702 V_mean = numpy.zeros(self.Num_Hei)
703 del_V = numpy.zeros(self.Num_Hei)
704 Z = numpy.zeros(self.Num_Hei)
705 Ze = numpy.zeros(self.Num_Hei)
706 RR = numpy.zeros(self.Num_Hei)
707
708 Range = dataOut.heightList*1000.
709
710 for R in range(self.Num_Hei):
711
712 h = Range[R] + Altitude #Range from ground to radar pulse altitude
713 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
714
715 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
716
717 '''NOTA: ETA(n) dn = ETA(f) df
718
719 dn = 1 Diferencial de muestreo
720 df = ETA(n) / ETA(f)
721
722 '''
723
724 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
725
726 ETAv[:,R]=ETAn[:,R]/dV
727
728 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
729
730 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
731
732 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
733
734 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
735
736 try:
737 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
738 except:
739 popt01=numpy.zeros(3)
740 popt01[1]= DMoments[1]
741
742 if popt01[1]<0 or popt01[1]>20:
743 popt01[1]=numpy.NaN
744
745
746 V_mean[R]=popt01[1]
747
748 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
749
750 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
751
752 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
753
754
755
756 RR2 = (Z/200)**(1/1.6)
757 dBRR = 10*numpy.log10(RR)
758 dBRR2 = 10*numpy.log10(RR2)
759
760 dBZe = 10*numpy.log10(Ze)
761 dBZ = 10*numpy.log10(Z)
743 ExpConstant = 10**(40/10) #Constante Experimental
744
745 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
746 for i in range(self.Num_Chn):
747 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
748 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
749
750 SPCmean = numpy.mean(SignalPower, 0)
751 Pr = SPCmean[:,:]/dataOut.normFactor
752
753 # Declaring auxiliary variables
754 Range = dataOut.heightList*1000. #Range in m
755 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
756 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
757 zMtrx = rMtrx+Altitude
758 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
759 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
760
761 # height dependence to air density Foote and Du Toit (1969)
762 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
763 VMtrx = VelMtrx / delv_z #Normalized velocity
764 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
765 # Diameter is related to the fall speed of falling drops
766 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
767 # Only valid for D>= 0.16 mm
768 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
769
770 #Calculate Radar Reflectivity ETAn
771 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
772 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
773 # Radar Cross Section
774 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
775 # Drop Size Distribution
776 DSD = ETAn / sigmaD
777 # Equivalente Reflectivy
778 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
779 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
780 # RainFall Rate
781 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
782
783 # Censoring the data
784 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
785 SNRth = 10**(SNRdBlimit/10) #-30dB
786 novalid = numpy.where((dataOut.data_snr[0,:] <SNRth) | (dataOut.data_snr[1,:] <SNRth) | (dataOut.data_snr[2,:] <SNRth)) # AND condition. Maybe OR condition better
787 W = numpy.nanmean(dataOut.data_dop,0)
788 W[novalid] = numpy.NaN
789 Ze_org[novalid] = numpy.NaN
790 RR[novalid] = numpy.NaN
762 791
763 792 dataOut.data_output = RR[8]
764 793 dataOut.data_param = numpy.ones([3,self.Num_Hei])
765 794 dataOut.channelList = [0,1,2]
766
767 dataOut.data_param[0]=dBZ
768 dataOut.data_param[1]=V_mean
795
796 dataOut.data_param[0]=10*numpy.log10(Ze_org)
797 dataOut.data_param[1]=-W
769 798 dataOut.data_param[2]=RR
770 799
800 # print ('Leaving PrecepitationProc ... ')
771 801 return dataOut
772 802
773 803 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
774 804
775 805 NPW = dataOut.NPW
776 806 COFA = dataOut.COFA
777 807
778 808 SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
779 809 RadarConst = dataOut.RadarConst
780 810 #frequency = 34.85*10**9
781 811
782 812 ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
783 813 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
784 814
785 815 ETA = numpy.sum(SNR,1)
786 816
787 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
817 ETA = numpy.where(ETA != 0. , ETA, numpy.NaN)
788 818
789 819 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
790 820
791 821 for r in range(self.Num_Hei):
792 822
793 823 Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
794 824 #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
795 825
796 826 return Ze
797 827
798 828 # def GetRadarConstant(self):
799 829 #
800 830 # """
801 831 # Constants:
802 832 #
803 833 # Pt: Transmission Power dB 5kW 5000
804 834 # Gt: Transmission Gain dB 24.7 dB 295.1209
805 835 # Gr: Reception Gain dB 18.5 dB 70.7945
806 836 # Lambda: Wavelenght m 0.6741 m 0.6741
807 837 # aL: Attenuation loses dB 4dB 2.5118
808 838 # tauW: Width of transmission pulse s 4us 4e-6
809 839 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
810 840 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
811 841 #
812 842 # """
813 843 #
814 844 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
815 845 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
816 846 # RadarConstant = Numerator / Denominator
817 847 #
818 848 # return RadarConstant
819 849
820 850
821 851
822 852 class FullSpectralAnalysis(Operation):
823 853
824 854 """
825 855 Function that implements Full Spectral Analysis technique.
826 856
827 857 Input:
828 858 self.dataOut.data_pre : SelfSpectra and CrossSpectra data
829 859 self.dataOut.groupList : Pairlist of channels
830 860 self.dataOut.ChanDist : Physical distance between receivers
831 861
832 862
833 863 Output:
834 864
835 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
865 self.dataOut.data_output : Zonal wind, Meridional wind, and Vertical wind
836 866
837 867
838 868 Parameters affected: Winds, height range, SNR
839 869
840 870 """
841 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7, minheight=None, maxheight=None):
842
843 self.indice=int(numpy.random.rand()*1000)
871 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
872 minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
844 873
845 874 spc = dataOut.data_pre[0].copy()
846 875 cspc = dataOut.data_pre[1]
847
848 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
849
850 SNRspc = spc.copy()
851 SNRspc[:,:,0:7]= numpy.NaN
852
853 """##########################################"""
854
855
856 nChannel = spc.shape[0]
857 nProfiles = spc.shape[1]
858 876 nHeights = spc.shape[2]
859 877
860 878 # first_height = 0.75 #km (ref: data header 20170822)
861 879 # resolution_height = 0.075 #km
862 880 '''
863 881 finding height range. check this when radar parameters are changed!
864 882 '''
865 883 if maxheight is not None:
866 884 # range_max = math.ceil((maxheight - first_height) / resolution_height) # theoretical
867 885 range_max = math.ceil(13.26 * maxheight - 3) # empirical, works better
868 886 else:
869 887 range_max = nHeights
870 888 if minheight is not None:
871 889 # range_min = int((minheight - first_height) / resolution_height) # theoretical
872 890 range_min = int(13.26 * minheight - 5) # empirical, works better
873 891 if range_min < 0:
874 892 range_min = 0
875 893 else:
876 894 range_min = 0
877 895
878 896 pairsList = dataOut.groupList
879 897 if dataOut.ChanDist is not None :
880 898 ChanDist = dataOut.ChanDist
881 899 else:
882 900 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
883 901
884 FrecRange = dataOut.spc_range[0]
885
886 data_SNR=numpy.zeros([nProfiles])
887 noise = dataOut.noise
888
889 dataOut.data_snr = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
890
891 dataOut.data_snr[numpy.where( dataOut.data_snr <0 )] = 1e-20
892
902 # 4 variables: zonal, meridional, vertical, and average SNR
903 data_param = numpy.zeros([4,nHeights]) * numpy.NaN
904 velocityX = numpy.zeros([nHeights]) * numpy.NaN
905 velocityY = numpy.zeros([nHeights]) * numpy.NaN
906 velocityZ = numpy.zeros([nHeights]) * numpy.NaN
893 907
894 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
895
896 velocityX=[]
897 velocityY=[]
898 velocityV=[]
899
900 dbSNR = 10*numpy.log10(dataOut.data_snr)
901 dbSNR = numpy.average(dbSNR,0)
908 dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
902 909
903 910 '''***********************************************WIND ESTIMATION**************************************'''
904
905 911 for Height in range(nHeights):
906 912
907 913 if Height >= range_min and Height < range_max:
908 # error_code unused, yet maybe useful for future analysis.
909 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
910 else:
911 Vzon,Vmer,Vver = 0., 0., numpy.NaN
912
913
914 if abs(Vzon) < 100. and abs(Vzon) > 0. and abs(Vmer) < 100. and abs(Vmer) > 0.:
915 velocityX=numpy.append(velocityX, Vzon)
916 velocityY=numpy.append(velocityY, -Vmer)
917
918 else:
919 velocityX=numpy.append(velocityX, numpy.NaN)
920 velocityY=numpy.append(velocityY, numpy.NaN)
921
922 if dbSNR[Height] > SNRlimit:
923 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version)
924 else:
925 velocityV=numpy.append(velocityV, numpy.NaN)
926
927
928 '''Change the numpy.array (velocityX) sign when trying to process BLTR data (Erick)'''
929 data_output[0] = numpy.array(velocityX)
930 data_output[1] = numpy.array(velocityY)
931 data_output[2] = velocityV
932
933
934 dataOut.data_output = data_output
935
914 # error_code will be useful in future analysis
915 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
916 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
917
918 if abs(Vzon) < 100. and abs(Vmer) < 100.:
919 velocityX[Height] = Vzon
920 velocityY[Height] = -Vmer
921 velocityZ[Height] = Vver
922
923 # Censoring data with SNR threshold
924 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
925
926 data_param[0] = velocityX
927 data_param[1] = velocityY
928 data_param[2] = velocityZ
929 data_param[3] = dbSNR
930 dataOut.data_param = data_param
936 931 return dataOut
937 932
938
939 933 def moving_average(self,x, N=2):
940 934 """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
941 935 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
942 936
943 937 def gaus(self,xSamples,Amp,Mu,Sigma):
944 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
938 return Amp * numpy.exp(-0.5*((xSamples - Mu)/Sigma)**2)
945 939
946 940 def Moments(self, ySamples, xSamples):
947 '''***
948 Variables corresponding to moments of distribution.
949 Also used as initial coefficients for curve_fit.
950 Vr was corrected. Only a velocity when x is velocity, of course.
951 ***'''
952 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
953 yNorm = ySamples / Pot
954 x_range = (numpy.max(xSamples)-numpy.min(xSamples))
955 Vr = numpy.nansum( yNorm * xSamples )*x_range/len(xSamples) # Velocidad radial, mu, corrimiento doppler, primer momento
956 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
957 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
958
959 return numpy.array([Pot, Vr, Desv])
941 Power = numpy.nanmean(ySamples) # Power, 0th Moment
942 yNorm = ySamples / numpy.nansum(ySamples)
943 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
944 Sigma2 = numpy.nansum(yNorm * (xSamples - RadVel)**2) # Spectral Width, 2nd Moment
945 StdDev = numpy.sqrt(numpy.abs(Sigma2)) # Desv. Estandar, Ancho espectral
946 return numpy.array([Power,RadVel,StdDev])
960 947
961 948 def StopWindEstimation(self, error_code):
962 '''
963 the wind calculation and returns zeros
964 '''
965 Vzon = 0
966 Vmer = 0
967 Vver = numpy.nan
949 Vzon = numpy.NaN
950 Vmer = numpy.NaN
951 Vver = numpy.NaN
968 952 return Vzon, Vmer, Vver, error_code
969 953
970 954 def AntiAliasing(self, interval, maxstep):
971 955 """
972 956 function to prevent errors from aliased values when computing phaseslope
973 957 """
974 antialiased = numpy.zeros(len(interval))*0.0
958 antialiased = numpy.zeros(len(interval))
975 959 copyinterval = interval.copy()
976 960
977 961 antialiased[0] = copyinterval[0]
978 962
979 963 for i in range(1,len(antialiased)):
980
981 964 step = interval[i] - interval[i-1]
982
983 965 if step > maxstep:
984 966 copyinterval -= 2*numpy.pi
985 967 antialiased[i] = copyinterval[i]
986
987 968 elif step < maxstep*(-1):
988 969 copyinterval += 2*numpy.pi
989 970 antialiased[i] = copyinterval[i]
990
991 971 else:
992 972 antialiased[i] = copyinterval[i].copy()
993 973
994 974 return antialiased
995 975
996 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
976 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit, NegativeLimit, PositiveLimit, radfreq):
997 977 """
998 978 Function that Calculates Zonal, Meridional and Vertical wind velocities.
999 979 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
1000 980
1001 981 Input:
1002 982 spc, cspc : self spectra and cross spectra data. In Briggs notation something like S_i*(S_i)_conj, (S_j)_conj respectively.
1003 983 pairsList : Pairlist of channels
1004 984 ChanDist : array of xi_ij and eta_ij
1005 985 Height : height at which data is processed
1006 986 noise : noise in [channels] format for specific height
1007 987 Abbsisarange : range of the frequencies or velocities
1008 988 dbSNR, SNRlimit : signal to noise ratio in db, lower limit
1009 989
1010 990 Output:
1011 991 Vzon, Vmer, Vver : wind velocities
1012 992 error_code : int that states where code is terminated
1013 993
1014 994 0 : no error detected
1015 995 1 : Gaussian of mean spc exceeds widthlimit
1016 996 2 : no Gaussian of mean spc found
1017 997 3 : SNR to low or velocity to high -> prec. e.g.
1018 998 4 : at least one Gaussian of cspc exceeds widthlimit
1019 999 5 : zero out of three cspc Gaussian fits converged
1020 1000 6 : phase slope fit could not be found
1021 1001 7 : arrays used to fit phase have different length
1022 1002 8 : frequency range is either too short (len <= 5) or very long (> 30% of cspc)
1023 1003
1024 1004 """
1025 1005
1026 1006 error_code = 0
1027 1007
1028
1029 SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height
1030 phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels
1031 CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values
1032 PhaseSlope = numpy.zeros(spc.shape[0]) # slope of the phases, channelwise
1033 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
1034 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1035 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
1036 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0]
1037
1038 SPCmoments_vel = self.Moments(SPCav, xVel ) # SPCmoments_vel[1] corresponds to vertical velocity and is used to determine if signal corresponds to wind (if .. <3)
1039 CSPCmoments = []
1040
1008 nChan = spc.shape[0]
1009 nProf = spc.shape[1]
1010 nPair = cspc.shape[0]
1011
1012 SPC_Samples = numpy.zeros([nChan, nProf]) # for normalized spc values for one height
1013 CSPC_Samples = numpy.zeros([nPair, nProf], dtype=numpy.complex_) # for normalized cspc values
1014 phase = numpy.zeros([nPair, nProf]) # phase between channels
1015 PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
1016 PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
1017 xFrec = AbbsisaRange[0][:-1] # frequency range
1018 xVel = AbbsisaRange[2][:-1] # velocity range
1019 xSamples = xFrec # the frequency range is taken
1020 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
1021
1022 # only consider velocities with in NegativeLimit and PositiveLimit
1023 if (NegativeLimit is None):
1024 NegativeLimit = numpy.min(xVel)
1025 if (PositiveLimit is None):
1026 PositiveLimit = numpy.max(xVel)
1027 xvalid = numpy.where((xVel > NegativeLimit) & (xVel < PositiveLimit))
1028 xSamples_zoom = xSamples[xvalid]
1041 1029
1042 1030 '''Getting Eij and Nij'''
1043
1044 1031 Xi01, Xi02, Xi12 = ChanDist[:,0]
1045 1032 Eta01, Eta02, Eta12 = ChanDist[:,1]
1046 1033
1047 # update nov 19
1048 widthlimit = 7 # maximum width in Hz of the gaussian, empirically determined. Anything above 10 is unrealistic, often values between 1 and 5 correspond to proper fits.
1049
1034 # spwd limit - updated by D. ScipiΓ³n 30.03.2021
1035 widthlimit = 10
1050 1036 '''************************* SPC is normalized ********************************'''
1051
1052 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1053 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1054
1055 for i in range(spc.shape[0]):
1056
1057 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1058
1059 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1060 normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm)
1061
1062 xSamples = xFrec # the frequency range is taken
1063 SPC_Samples[i] = normalized_spc # Normalized SPC values are taken
1037 spc_norm = spc.copy()
1038 # For each channel
1039 for i in range(nChan):
1040 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
1041 SPC_Samples[i] = spc_sub / (numpy.nansum(spc_sub) * delta_x)
1064 1042
1065 1043 '''********************** FITTING MEAN SPC GAUSSIAN **********************'''
1066 1044
1067 1045 """ the gaussian of the mean: first subtract noise, then normalize. this is legal because
1068 1046 you only fit the curve and don't need the absolute value of height for calculation,
1069 1047 only for estimation of width. for normalization of cross spectra, you need initial,
1070 1048 unnormalized self-spectra With noise.
1071 1049
1072 1050 Technically, you don't even need to normalize the self-spectra, as you only need the
1073 1051 width of the peak. However, it was left this way. Note that the normalization has a flaw:
1074 1052 due to subtraction of the noise, some values are below zero. Raw "spc" values should be
1075 1053 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1076 1054 """
1077
1078 SPCMean = numpy.average(SPC_Samples, axis=0)
1079
1080 popt = [1e-10,0,1e-10]
1081 SPCMoments = self.Moments(SPCMean, xSamples)
1082
1083 if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3:
1055 # initial conditions
1056 popt = [1e-10,0,1e-10]
1057 # Spectra average
1058 SPCMean = numpy.average(SPC_Samples,0)
1059 # Moments in frequency
1060 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1061
1062 # Gauss Fit SPC in frequency domain
1063 if dbSNR > SNRlimit: # only if SNR > SNRth
1084 1064 try:
1085 popt,pcov = curve_fit(self.gaus,xSamples,SPCMean,p0=SPCMoments)#, bounds=(-numpy.inf, [numpy.inf, numpy.inf, 10])). Setting bounds does not make the code faster but only keeps the fit from finding the minimum.
1086 if popt[2] > widthlimit: # CONDITION
1065 popt,pcov = curve_fit(self.gaus,xSamples_zoom,SPCMean[xvalid],p0=SPCMoments)
1066 if popt[2] <= 0 or popt[2] > widthlimit: # CONDITION
1087 1067 return self.StopWindEstimation(error_code = 1)
1088
1089 FitGauss = self.gaus(xSamples,*popt)
1090
1068 FitGauss = self.gaus(xSamples_zoom,*popt)
1091 1069 except :#RuntimeError:
1092 1070 return self.StopWindEstimation(error_code = 2)
1093
1094 1071 else:
1095 1072 return self.StopWindEstimation(error_code = 3)
1096 1073
1097
1098
1099 1074 '''***************************** CSPC Normalization *************************
1100 new section:
1101 1075 The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
1102 1076 influence the norm which is not desired. First, a range is identified where the
1103 1077 wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
1104 1078 around it gets cut off and values replaced by mean determined by the boundary
1105 1079 data -> sum_noise (spc is not normalized here, thats why the noise is important)
1106 1080
1107 1081 The sums are then added and multiplied by range/datapoints, because you need
1108 1082 an integral and not a sum for normalization.
1109 1083
1110 1084 A norm is found according to Briggs 92.
1111 1085 '''
1112
1113 radarWavelength = 0.6741 # meters
1114 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1115 # count_limit_freq = numpy.max(xFrec)
1116
1117 channel_integrals = numpy.zeros(3)
1118
1119 for i in range(spc.shape[0]):
1120 '''
1121 find the point in array corresponding to count_limit frequency.
1122 sum over all frequencies in the range around zero Hz @ math.ceil(N_freq/2)
1123 '''
1124 N_freq = numpy.count_nonzero(~numpy.isnan(spc[i,:]))
1125 count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point
1126 sum_wind = numpy.nansum( spc[i, (math.ceil(N_freq/2) - count_limit_int) : (math.ceil(N_freq / 2) + count_limit_int)] ) #N_freq/2 is where frequency (velocity) is zero, i.e. middle of spectrum.
1127 sum_noise = (numpy.mean(spc[i, :4]) + numpy.mean(spc[i, -6:-2]))/2.0 * (N_freq - 2*count_limit_int)
1128 channel_integrals[i] = (sum_noise + sum_wind) * (2*numpy.max(xFrec) / N_freq)
1129
1130
1131 cross_integrals_peak = numpy.zeros(3)
1132 # cross_integrals_totalrange = numpy.zeros(3)
1133
1134 for i in range(spc.shape[0]):
1135
1136 cspc_norm = cspc[i,:].copy() # cspc not smoothed here or in previous version
1137
1086 # for each pair
1087 for i in range(nPair):
1088 cspc_norm = cspc[i,:].copy()
1138 1089 chan_index0 = pairsList[i][0]
1139 1090 chan_index1 = pairsList[i][1]
1140
1141 cross_integrals_peak[i] = channel_integrals[chan_index0]*channel_integrals[chan_index1]
1142 normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i])
1143 CSPC_Samples[i] = normalized_cspc
1144
1145 ''' Finding cross integrals without subtracting any peaks:'''
1146 # FactorNorm0 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index0,:]))
1147 # FactorNorm1 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index1,:]))
1148 # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1
1149 # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i])
1150 # CSPC_Samples[i] = normalized_cspc
1151
1152
1091 CSPC_Samples[i] = cspc_norm / (numpy.sqrt(numpy.nansum(spc_norm[chan_index0])*numpy.nansum(spc_norm[chan_index1])) * delta_x)
1153 1092 phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
1154 1093
1094 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0,xvalid]), xSamples_zoom),
1095 self.Moments(numpy.abs(CSPC_Samples[1,xvalid]), xSamples_zoom),
1096 self.Moments(numpy.abs(CSPC_Samples[2,xvalid]), xSamples_zoom)])
1155 1097
1156 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0]), xSamples),
1157 self.Moments(numpy.abs(CSPC_Samples[1]), xSamples),
1158 self.Moments(numpy.abs(CSPC_Samples[2]), xSamples)])
1159
1160
1161 '''***Sorting out NaN entries***'''
1162 CSPCMask01 = numpy.abs(CSPC_Samples[0])
1163 CSPCMask02 = numpy.abs(CSPC_Samples[1])
1164 CSPCMask12 = numpy.abs(CSPC_Samples[2])
1165
1166 mask01 = ~numpy.isnan(CSPCMask01)
1167 mask02 = ~numpy.isnan(CSPCMask02)
1168 mask12 = ~numpy.isnan(CSPCMask12)
1169
1170 CSPCMask01 = CSPCMask01[mask01]
1171 CSPCMask02 = CSPCMask02[mask02]
1172 CSPCMask12 = CSPCMask12[mask12]
1173
1174
1175 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1176 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1098 popt01, popt02, popt12 = [1e-10,0,1e-10], [1e-10,0,1e-10] ,[1e-10,0,1e-10]
1099 FitGauss01, FitGauss02, FitGauss12 = numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples))
1177 1100
1178 1101 '''*******************************FIT GAUSS CSPC************************************'''
1179
1180 1102 try:
1181 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1103 popt01,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[0][xvalid]),p0=CSPCmoments[0])
1182 1104 if popt01[2] > widthlimit: # CONDITION
1183 1105 return self.StopWindEstimation(error_code = 4)
1184
1185 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1106 popt02,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[1][xvalid]),p0=CSPCmoments[1])
1186 1107 if popt02[2] > widthlimit: # CONDITION
1187 1108 return self.StopWindEstimation(error_code = 4)
1188
1189 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1109 popt12,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[2][xvalid]),p0=CSPCmoments[2])
1190 1110 if popt12[2] > widthlimit: # CONDITION
1191 1111 return self.StopWindEstimation(error_code = 4)
1192 1112
1193 FitGauss01 = self.gaus(xSamples, *popt01)
1194 FitGauss02 = self.gaus(xSamples, *popt02)
1195 FitGauss12 = self.gaus(xSamples, *popt12)
1196
1113 FitGauss01 = self.gaus(xSamples_zoom, *popt01)
1114 FitGauss02 = self.gaus(xSamples_zoom, *popt02)
1115 FitGauss12 = self.gaus(xSamples_zoom, *popt12)
1197 1116 except:
1198 1117 return self.StopWindEstimation(error_code = 5)
1199 1118
1200 1119
1201 1120 '''************* Getting Fij ***************'''
1202
1203
1204 #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located
1205 # -> PointGauCenter
1121 # x-axis point of the gaussian where the center is located from GaussFit of spectra
1206 1122 GaussCenter = popt[1]
1207 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1208 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1123 ClosestCenter = xSamples_zoom[numpy.abs(xSamples_zoom-GaussCenter).argmin()]
1124 PointGauCenter = numpy.where(xSamples_zoom==ClosestCenter)[0][0]
1209 1125
1210 #Punto e^-1 hubicado en la Gaussiana -- point where e^-1 is located in the gaussian
1126 # Point where e^-1 is located in the gaussian
1211 1127 PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1)
1212 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1128 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # The closest point to"Peminus1" in "FitGauss"
1213 1129 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1214
1215 Fij = numpy.abs(xSamples[PointFij] - xSamples[PointGauCenter])
1130 Fij = numpy.abs(xSamples_zoom[PointFij] - xSamples_zoom[PointGauCenter])
1216 1131
1217 1132 '''********** Taking frequency ranges from mean SPCs **********'''
1218
1219 #GaussCenter = popt[1] #Primer momento 01
1220 GauWidth = popt[2] * 3/2 #Ancho de banda de Gau01 -- Bandwidth of Gau01 TODO why *3/2?
1133 GauWidth = popt[2] * 3/2 # Bandwidth of Gau01
1221 1134 Range = numpy.empty(2)
1222 1135 Range[0] = GaussCenter - GauWidth
1223 1136 Range[1] = GaussCenter + GauWidth
1224 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) -- Point in x-axis where the bandwidth is located (min:max)
1225 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1226 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1227
1228 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1229 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1230
1137 # Point in x-axis where the bandwidth is located (min:max)
1138 ClosRangeMin = xSamples_zoom[numpy.abs(xSamples_zoom-Range[0]).argmin()]
1139 ClosRangeMax = xSamples_zoom[numpy.abs(xSamples_zoom-Range[1]).argmin()]
1140 PointRangeMin = numpy.where(xSamples_zoom==ClosRangeMin)[0][0]
1141 PointRangeMax = numpy.where(xSamples_zoom==ClosRangeMax)[0][0]
1231 1142 Range = numpy.array([ PointRangeMin, PointRangeMax ])
1232
1233 FrecRange = xFrec[ Range[0] : Range[1] ]
1234
1143 FrecRange = xSamples_zoom[ Range[0] : Range[1] ]
1235 1144
1236 1145 '''************************** Getting Phase Slope ***************************'''
1237
1238 for i in range(1,3): # Changed to only compute two
1239
1240 if len(FrecRange) > 5 and len(FrecRange) < spc.shape[1] * 0.3:
1241 # PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=1) #used before to smooth phase with N=3
1242 PhaseRange = phase[i,Range[0]:Range[1]].copy()
1243
1146 for i in range(nPair):
1147 if len(FrecRange) > 5:
1148 PhaseRange = phase[i, xvalid[0][Range[0]:Range[1]]].copy()
1244 1149 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1245
1246
1247 1150 if len(FrecRange) == len(PhaseRange):
1248
1249 1151 try:
1250 1152 slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
1251 1153 PhaseSlope[i] = slope
1252 1154 PhaseInter[i] = intercept
1253
1254 1155 except:
1255 1156 return self.StopWindEstimation(error_code = 6)
1256
1257 1157 else:
1258 1158 return self.StopWindEstimation(error_code = 7)
1259
1260 1159 else:
1261 1160 return self.StopWindEstimation(error_code = 8)
1262 1161
1263
1264
1265 1162 '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
1266 1163
1267 1164 '''Getting constant C'''
1268 1165 cC=(Fij*numpy.pi)**2
1269 1166
1270 1167 '''****** Getting constants F and G ******'''
1271 1168 MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1272 MijResult0 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1273 MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1274 MijResults = numpy.array([MijResult0,MijResult1])
1169 # MijEijNij = numpy.array([[Xi01,Eta01], [Xi02,Eta02], [Xi12,Eta12]])
1170 # MijResult0 = (-PhaseSlope[0] * cC) / (2*numpy.pi)
1171 MijResult1 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1172 MijResult2 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1173 # MijResults = numpy.array([MijResult0, MijResult1, MijResult2])
1174 MijResults = numpy.array([MijResult1, MijResult2])
1275 1175 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1276 1176
1277 1177 '''****** Getting constants A, B and H ******'''
1278 1178 W01 = numpy.nanmax( FitGauss01 )
1279 1179 W02 = numpy.nanmax( FitGauss02 )
1280 1180 W12 = numpy.nanmax( FitGauss12 )
1281 1181
1282 WijResult0 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1283 WijResult1 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1284 WijResult2 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1285
1286 WijResults = numpy.array([WijResult0, WijResult1, WijResult2])
1182 WijResult01 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1183 WijResult02 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1184 WijResult12 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1185 WijResults = numpy.array([WijResult01, WijResult02, WijResult12])
1287 1186
1288 1187 WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1289 1188 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1290 1189
1291 1190 VxVy = numpy.array([[cA,cH],[cH,cB]])
1292 1191 VxVyResults = numpy.array([-cF,-cG])
1293 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1294
1295 Vzon = Vy
1296 Vmer = Vx
1297
1298 # Vmag=numpy.sqrt(Vzon**2+Vmer**2) # unused
1299 # Vang=numpy.arctan2(Vmer,Vzon) # unused
1300
1301
1302 ''' using frequency as abscissa. Due to three channels, the offzenith angle is zero
1303 and Vrad equal to Vver. formula taken from Briggs 92, figure 4.
1304 '''
1305 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange) > 4:
1306 Vver = 0.5 * radarWavelength * popt[1] * 100 # *100 to get cm (/s)
1307 else:
1308 Vver = numpy.NaN
1309
1192 (Vmer,Vzon) = numpy.linalg.solve(VxVy, VxVyResults)
1193 Vver = -SPCMoments[1]*SPEED_OF_LIGHT/(2*radfreq)
1310 1194 error_code = 0
1311 1195
1312 1196 return Vzon, Vmer, Vver, error_code
1313 1197
1314
1315 1198 class SpectralMoments(Operation):
1316 1199
1317 1200 '''
1318 1201 Function SpectralMoments()
1319 1202
1320 1203 Calculates moments (power, mean, standard deviation) and SNR of the signal
1321 1204
1322 1205 Type of dataIn: Spectra
1323 1206
1324 1207 Configuration Parameters:
1325 1208
1326 1209 dirCosx : Cosine director in X axis
1327 1210 dirCosy : Cosine director in Y axis
1328 1211
1329 1212 elevation :
1330 1213 azimuth :
1331 1214
1332 1215 Input:
1333 1216 channelList : simple channel list to select e.g. [2,3,7]
1334 1217 self.dataOut.data_pre : Spectral data
1335 1218 self.dataOut.abscissaList : List of frequencies
1336 1219 self.dataOut.noise : Noise level per channel
1337 1220
1338 1221 Affected:
1339 1222 self.dataOut.moments : Parameters per channel
1340 1223 self.dataOut.data_snr : SNR per channel
1341 1224
1342 1225 '''
1343 1226
1344 1227 def run(self, dataOut):
1345 1228
1346 1229 data = dataOut.data_pre[0]
1347 1230 absc = dataOut.abscissaList[:-1]
1348 1231 noise = dataOut.noise
1349 1232 nChannel = data.shape[0]
1350 1233 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
1351 1234
1352 1235 for ind in range(nChannel):
1353 1236 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1354 1237
1355 1238 dataOut.moments = data_param[:,1:,:]
1356 1239 dataOut.data_snr = data_param[:,0]
1357 1240 dataOut.data_pow = data_param[:,1]
1358 1241 dataOut.data_dop = data_param[:,2]
1359 1242 dataOut.data_width = data_param[:,3]
1360 1243
1361 1244 return dataOut
1362 1245
1363 1246 def __calculateMoments(self, oldspec, oldfreq, n0,
1364 1247 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
1365 1248
1366 1249 if (nicoh is None): nicoh = 1
1367 1250 if (graph is None): graph = 0
1368 1251 if (smooth is None): smooth = 0
1369 1252 elif (self.smooth < 3): smooth = 0
1370 1253
1371 1254 if (type1 is None): type1 = 0
1372 1255 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
1373 1256 if (snrth is None): snrth = -3
1374 1257 if (dc is None): dc = 0
1375 1258 if (aliasing is None): aliasing = 0
1376 1259 if (oldfd is None): oldfd = 0
1377 1260 if (wwauto is None): wwauto = 0
1378 1261
1379 1262 if (n0 < 1.e-20): n0 = 1.e-20
1380 1263
1381 1264 freq = oldfreq
1382 1265 vec_power = numpy.zeros(oldspec.shape[1])
1383 1266 vec_fd = numpy.zeros(oldspec.shape[1])
1384 1267 vec_w = numpy.zeros(oldspec.shape[1])
1385 1268 vec_snr = numpy.zeros(oldspec.shape[1])
1386 1269
1387 1270 # oldspec = numpy.ma.masked_invalid(oldspec)
1388 1271
1389 1272 for ind in range(oldspec.shape[1]):
1390 1273
1391 1274 spec = oldspec[:,ind]
1392 1275 aux = spec*fwindow
1393 1276 max_spec = aux.max()
1394 1277 m = aux.tolist().index(max_spec)
1395 1278
1396 #Smooth
1279 # Smooth
1397 1280 if (smooth == 0):
1398 1281 spec2 = spec
1399 1282 else:
1400 1283 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1401 1284
1402 # Calculo de Momentos
1285 # Moments Estimation
1403 1286 bb = spec2[numpy.arange(m,spec2.size)]
1404 1287 bb = (bb<n0).nonzero()
1405 1288 bb = bb[0]
1406 1289
1407 1290 ss = spec2[numpy.arange(0,m + 1)]
1408 1291 ss = (ss<n0).nonzero()
1409 1292 ss = ss[0]
1410 1293
1411 1294 if (bb.size == 0):
1412 1295 bb0 = spec.size - 1 - m
1413 1296 else:
1414 1297 bb0 = bb[0] - 1
1415 1298 if (bb0 < 0):
1416 1299 bb0 = 0
1417 1300
1418 1301 if (ss.size == 0):
1419 1302 ss1 = 1
1420 1303 else:
1421 1304 ss1 = max(ss) + 1
1422 1305
1423 1306 if (ss1 > m):
1424 1307 ss1 = m
1425 1308
1426 1309 valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
1427 1310
1428 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1311 signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. ScipiΓ³n added with correct definition
1312 total_power = (spec2[valid] * fwindow[valid]).mean() # D. ScipiΓ³n added with correct definition
1313 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1429 1314 fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
1430 1315 w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
1431 1316 snr = (spec2.mean()-n0)/n0
1432 1317 if (snr < 1.e-20) :
1433 1318 snr = 1.e-20
1434 1319
1435 vec_power[ind] = power
1320 # vec_power[ind] = power #D. ScipiΓ³n replaced with the line below
1321 vec_power[ind] = total_power
1436 1322 vec_fd[ind] = fd
1437 1323 vec_w[ind] = w
1438 1324 vec_snr[ind] = snr
1439 1325
1440 1326 return numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1441 1327
1442 1328 #------------------ Get SA Parameters --------------------------
1443 1329
1444 1330 def GetSAParameters(self):
1445 1331 #SA en frecuencia
1446 1332 pairslist = self.dataOut.groupList
1447 1333 num_pairs = len(pairslist)
1448 1334
1449 1335 vel = self.dataOut.abscissaList
1450 1336 spectra = self.dataOut.data_pre
1451 1337 cspectra = self.dataIn.data_cspc
1452 1338 delta_v = vel[1] - vel[0]
1453 1339
1454 1340 #Calculating the power spectrum
1455 1341 spc_pow = numpy.sum(spectra, 3)*delta_v
1456 1342 #Normalizing Spectra
1457 1343 norm_spectra = spectra/spc_pow
1458 1344 #Calculating the norm_spectra at peak
1459 1345 max_spectra = numpy.max(norm_spectra, 3)
1460 1346
1461 1347 #Normalizing Cross Spectra
1462 1348 norm_cspectra = numpy.zeros(cspectra.shape)
1463 1349
1464 1350 for i in range(num_chan):
1465 1351 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
1466 1352
1467 1353 max_cspectra = numpy.max(norm_cspectra,2)
1468 1354 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
1469 1355
1470 1356 for i in range(num_pairs):
1471 1357 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
1472 1358 #------------------- Get Lags ----------------------------------
1473 1359
1474 1360 class SALags(Operation):
1475 1361 '''
1476 1362 Function GetMoments()
1477 1363
1478 1364 Input:
1479 1365 self.dataOut.data_pre
1480 1366 self.dataOut.abscissaList
1481 1367 self.dataOut.noise
1482 1368 self.dataOut.normFactor
1483 1369 self.dataOut.data_snr
1484 1370 self.dataOut.groupList
1485 1371 self.dataOut.nChannels
1486 1372
1487 1373 Affected:
1488 1374 self.dataOut.data_param
1489 1375
1490 1376 '''
1491 1377 def run(self, dataOut):
1492 1378 data_acf = dataOut.data_pre[0]
1493 1379 data_ccf = dataOut.data_pre[1]
1494 1380 normFactor_acf = dataOut.normFactor[0]
1495 1381 normFactor_ccf = dataOut.normFactor[1]
1496 1382 pairs_acf = dataOut.groupList[0]
1497 1383 pairs_ccf = dataOut.groupList[1]
1498 1384
1499 1385 nHeights = dataOut.nHeights
1500 1386 absc = dataOut.abscissaList
1501 1387 noise = dataOut.noise
1502 1388 SNR = dataOut.data_snr
1503 1389 nChannels = dataOut.nChannels
1504 1390 # pairsList = dataOut.groupList
1505 1391 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1506 1392
1507 1393 for l in range(len(pairs_acf)):
1508 1394 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
1509 1395
1510 1396 for l in range(len(pairs_ccf)):
1511 1397 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
1512 1398
1513 1399 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
1514 1400 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
1515 1401 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
1516 1402 return
1517 1403
1518 1404 # def __getPairsAutoCorr(self, pairsList, nChannels):
1519 1405 #
1520 1406 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1521 1407 #
1522 1408 # for l in range(len(pairsList)):
1523 1409 # firstChannel = pairsList[l][0]
1524 1410 # secondChannel = pairsList[l][1]
1525 1411 #
1526 1412 # #Obteniendo pares de Autocorrelacion
1527 1413 # if firstChannel == secondChannel:
1528 1414 # pairsAutoCorr[firstChannel] = int(l)
1529 1415 #
1530 1416 # pairsAutoCorr = pairsAutoCorr.astype(int)
1531 1417 #
1532 1418 # pairsCrossCorr = range(len(pairsList))
1533 1419 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1534 1420 #
1535 1421 # return pairsAutoCorr, pairsCrossCorr
1536 1422
1537 1423 def __calculateTaus(self, data_acf, data_ccf, lagRange):
1538 1424
1539 1425 lag0 = data_acf.shape[1]/2
1540 1426 #Funcion de Autocorrelacion
1541 1427 mean_acf = stats.nanmean(data_acf, axis = 0)
1542 1428
1543 1429 #Obtencion Indice de TauCross
1544 1430 ind_ccf = data_ccf.argmax(axis = 1)
1545 1431 #Obtencion Indice de TauAuto
1546 1432 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
1547 1433 ccf_lag0 = data_ccf[:,lag0,:]
1548 1434
1549 1435 for i in range(ccf_lag0.shape[0]):
1550 1436 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
1551 1437
1552 1438 #Obtencion de TauCross y TauAuto
1553 1439 tau_ccf = lagRange[ind_ccf]
1554 1440 tau_acf = lagRange[ind_acf]
1555 1441
1556 1442 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
1557 1443
1558 1444 tau_ccf[Nan1,Nan2] = numpy.nan
1559 1445 tau_acf[Nan1,Nan2] = numpy.nan
1560 1446 tau = numpy.vstack((tau_ccf,tau_acf))
1561 1447
1562 1448 return tau
1563 1449
1564 1450 def __calculateLag1Phase(self, data, lagTRange):
1565 1451 data1 = stats.nanmean(data, axis = 0)
1566 1452 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
1567 1453
1568 1454 phase = numpy.angle(data1[lag1,:])
1569 1455
1570 1456 return phase
1571 1457
1572 1458 class SpectralFitting(Operation):
1573 1459 '''
1574 1460 Function GetMoments()
1575 1461
1576 1462 Input:
1577 1463 Output:
1578 1464 Variables modified:
1579 1465 '''
1580 1466
1581 1467 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1582 1468
1583 1469
1584 1470 if path != None:
1585 1471 sys.path.append(path)
1586 1472 self.dataOut.library = importlib.import_module(file)
1587 1473
1588 1474 #To be inserted as a parameter
1589 1475 groupArray = numpy.array(groupList)
1590 1476 # groupArray = numpy.array([[0,1],[2,3]])
1591 1477 self.dataOut.groupList = groupArray
1592 1478
1593 1479 nGroups = groupArray.shape[0]
1594 1480 nChannels = self.dataIn.nChannels
1595 1481 nHeights=self.dataIn.heightList.size
1596 1482
1597 1483 #Parameters Array
1598 1484 self.dataOut.data_param = None
1599 1485
1600 1486 #Set constants
1601 1487 constants = self.dataOut.library.setConstants(self.dataIn)
1602 1488 self.dataOut.constants = constants
1603 1489 M = self.dataIn.normFactor
1604 1490 N = self.dataIn.nFFTPoints
1605 1491 ippSeconds = self.dataIn.ippSeconds
1606 1492 K = self.dataIn.nIncohInt
1607 1493 pairsArray = numpy.array(self.dataIn.pairsList)
1608 1494
1609 1495 #List of possible combinations
1610 1496 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1611 1497 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1612 1498
1613 1499 if getSNR:
1614 1500 listChannels = groupArray.reshape((groupArray.size))
1615 1501 listChannels.sort()
1616 1502 noise = self.dataIn.getNoise()
1617 1503 self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1618 1504
1619 1505 for i in range(nGroups):
1620 1506 coord = groupArray[i,:]
1621 1507
1622 1508 #Input data array
1623 1509 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1624 1510 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1625 1511
1626 1512 #Cross Spectra data array for Covariance Matrixes
1627 1513 ind = 0
1628 1514 for pairs in listComb:
1629 1515 pairsSel = numpy.array([coord[x],coord[y]])
1630 1516 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1631 1517 ind += 1
1632 1518 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1633 1519 dataCross = dataCross**2/K
1634 1520
1635 1521 for h in range(nHeights):
1636 1522
1637 1523 #Input
1638 1524 d = data[:,h]
1639 1525
1640 1526 #Covariance Matrix
1641 1527 D = numpy.diag(d**2/K)
1642 1528 ind = 0
1643 1529 for pairs in listComb:
1644 1530 #Coordinates in Covariance Matrix
1645 1531 x = pairs[0]
1646 1532 y = pairs[1]
1647 1533 #Channel Index
1648 1534 S12 = dataCross[ind,:,h]
1649 1535 D12 = numpy.diag(S12)
1650 1536 #Completing Covariance Matrix with Cross Spectras
1651 1537 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1652 1538 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1653 1539 ind += 1
1654 1540 Dinv=numpy.linalg.inv(D)
1655 1541 L=numpy.linalg.cholesky(Dinv)
1656 1542 LT=L.T
1657 1543
1658 1544 dp = numpy.dot(LT,d)
1659 1545
1660 1546 #Initial values
1661 1547 data_spc = self.dataIn.data_spc[coord,:,h]
1662 1548
1663 1549 if (h>0)and(error1[3]<5):
1664 1550 p0 = self.dataOut.data_param[i,:,h-1]
1665 1551 else:
1666 1552 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1667 1553
1668 1554 try:
1669 1555 #Least Squares
1670 1556 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1671 1557 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1672 1558 #Chi square error
1673 1559 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1674 1560 #Error with Jacobian
1675 1561 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1676 1562 except:
1677 1563 minp = p0*numpy.nan
1678 1564 error0 = numpy.nan
1679 1565 error1 = p0*numpy.nan
1680 1566
1681 1567 #Save
1682 1568 if self.dataOut.data_param is None:
1683 1569 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1684 1570 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1685 1571
1686 1572 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1687 1573 self.dataOut.data_param[i,:,h] = minp
1688 1574 return
1689 1575
1690 1576 def __residFunction(self, p, dp, LT, constants):
1691 1577
1692 1578 fm = self.dataOut.library.modelFunction(p, constants)
1693 1579 fmp=numpy.dot(LT,fm)
1694 1580
1695 1581 return dp-fmp
1696 1582
1697 1583 def __getSNR(self, z, noise):
1698 1584
1699 1585 avg = numpy.average(z, axis=1)
1700 1586 SNR = (avg.T-noise)/noise
1701 1587 SNR = SNR.T
1702 1588 return SNR
1703 1589
1704 1590 def __chisq(p,chindex,hindex):
1705 1591 #similar to Resid but calculates CHI**2
1706 1592 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1707 1593 dp=numpy.dot(LT,d)
1708 1594 fmp=numpy.dot(LT,fm)
1709 1595 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1710 1596 return chisq
1711 1597
1712 1598 class WindProfiler(Operation):
1713 1599
1714 1600 __isConfig = False
1715 1601
1716 1602 __initime = None
1717 1603 __lastdatatime = None
1718 1604 __integrationtime = None
1719 1605
1720 1606 __buffer = None
1721 1607
1722 1608 __dataReady = False
1723 1609
1724 1610 __firstdata = None
1725 1611
1726 1612 n = None
1727 1613
1728 1614 def __init__(self):
1729 1615 Operation.__init__(self)
1730 1616
1731 1617 def __calculateCosDir(self, elev, azim):
1732 1618 zen = (90 - elev)*numpy.pi/180
1733 1619 azim = azim*numpy.pi/180
1734 1620 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1735 1621 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1736 1622
1737 1623 signX = numpy.sign(numpy.cos(azim))
1738 1624 signY = numpy.sign(numpy.sin(azim))
1739 1625
1740 1626 cosDirX = numpy.copysign(cosDirX, signX)
1741 1627 cosDirY = numpy.copysign(cosDirY, signY)
1742 1628 return cosDirX, cosDirY
1743 1629
1744 1630 def __calculateAngles(self, theta_x, theta_y, azimuth):
1745 1631
1746 1632 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1747 1633 zenith_arr = numpy.arccos(dir_cosw)
1748 1634 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1749 1635
1750 1636 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1751 1637 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1752 1638
1753 1639 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1754 1640
1755 1641 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1756 1642
1757 1643 #
1758 1644 if horOnly:
1759 1645 A = numpy.c_[dir_cosu,dir_cosv]
1760 1646 else:
1761 1647 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1762 1648 A = numpy.asmatrix(A)
1763 1649 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1764 1650
1765 1651 return A1
1766 1652
1767 1653 def __correctValues(self, heiRang, phi, velRadial, SNR):
1768 1654 listPhi = phi.tolist()
1769 1655 maxid = listPhi.index(max(listPhi))
1770 1656 minid = listPhi.index(min(listPhi))
1771 1657
1772 1658 rango = list(range(len(phi)))
1773 1659 # rango = numpy.delete(rango,maxid)
1774 1660
1775 1661 heiRang1 = heiRang*math.cos(phi[maxid])
1776 1662 heiRangAux = heiRang*math.cos(phi[minid])
1777 1663 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1778 1664 heiRang1 = numpy.delete(heiRang1,indOut)
1779 1665
1780 1666 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1781 1667 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1782 1668
1783 1669 for i in rango:
1784 1670 x = heiRang*math.cos(phi[i])
1785 1671 y1 = velRadial[i,:]
1786 1672 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1787 1673
1788 1674 x1 = heiRang1
1789 1675 y11 = f1(x1)
1790 1676
1791 1677 y2 = SNR[i,:]
1792 1678 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1793 1679 y21 = f2(x1)
1794 1680
1795 1681 velRadial1[i,:] = y11
1796 1682 SNR1[i,:] = y21
1797 1683
1798 1684 return heiRang1, velRadial1, SNR1
1799 1685
1800 1686 def __calculateVelUVW(self, A, velRadial):
1801 1687
1802 1688 #Operacion Matricial
1803 1689 # velUVW = numpy.zeros((velRadial.shape[1],3))
1804 1690 # for ind in range(velRadial.shape[1]):
1805 1691 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1806 1692 # velUVW = velUVW.transpose()
1807 1693 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1808 1694 velUVW[:,:] = numpy.dot(A,velRadial)
1809 1695
1810 1696
1811 1697 return velUVW
1812 1698
1813 1699 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1814 1700
1815 1701 def techniqueDBS(self, kwargs):
1816 1702 """
1817 1703 Function that implements Doppler Beam Swinging (DBS) technique.
1818 1704
1819 1705 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1820 1706 Direction correction (if necessary), Ranges and SNR
1821 1707
1822 1708 Output: Winds estimation (Zonal, Meridional and Vertical)
1823 1709
1824 1710 Parameters affected: Winds, height range, SNR
1825 1711 """
1826 1712 velRadial0 = kwargs['velRadial']
1827 1713 heiRang = kwargs['heightList']
1828 1714 SNR0 = kwargs['SNR']
1829 1715
1830 1716 if 'dirCosx' in kwargs and 'dirCosy' in kwargs:
1831 1717 theta_x = numpy.array(kwargs['dirCosx'])
1832 1718 theta_y = numpy.array(kwargs['dirCosy'])
1833 1719 else:
1834 1720 elev = numpy.array(kwargs['elevation'])
1835 1721 azim = numpy.array(kwargs['azimuth'])
1836 1722 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1837 1723 azimuth = kwargs['correctAzimuth']
1838 1724 if 'horizontalOnly' in kwargs:
1839 1725 horizontalOnly = kwargs['horizontalOnly']
1840 1726 else: horizontalOnly = False
1841 1727 if 'correctFactor' in kwargs:
1842 1728 correctFactor = kwargs['correctFactor']
1843 1729 else: correctFactor = 1
1844 1730 if 'channelList' in kwargs:
1845 1731 channelList = kwargs['channelList']
1846 1732 if len(channelList) == 2:
1847 1733 horizontalOnly = True
1848 1734 arrayChannel = numpy.array(channelList)
1849 1735 param = param[arrayChannel,:,:]
1850 1736 theta_x = theta_x[arrayChannel]
1851 1737 theta_y = theta_y[arrayChannel]
1852 1738
1853 1739 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1854 1740 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
1855 1741 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1856 1742
1857 1743 #Calculo de Componentes de la velocidad con DBS
1858 1744 winds = self.__calculateVelUVW(A,velRadial1)
1859 1745
1860 1746 return winds, heiRang1, SNR1
1861 1747
1862 1748 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
1863 1749
1864 1750 nPairs = len(pairs_ccf)
1865 1751 posx = numpy.asarray(posx)
1866 1752 posy = numpy.asarray(posy)
1867 1753
1868 1754 #Rotacion Inversa para alinear con el azimuth
1869 1755 if azimuth!= None:
1870 1756 azimuth = azimuth*math.pi/180
1871 1757 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1872 1758 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1873 1759 else:
1874 1760 posx1 = posx
1875 1761 posy1 = posy
1876 1762
1877 1763 #Calculo de Distancias
1878 1764 distx = numpy.zeros(nPairs)
1879 1765 disty = numpy.zeros(nPairs)
1880 1766 dist = numpy.zeros(nPairs)
1881 1767 ang = numpy.zeros(nPairs)
1882 1768
1883 1769 for i in range(nPairs):
1884 1770 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
1885 1771 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
1886 1772 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1887 1773 ang[i] = numpy.arctan2(disty[i],distx[i])
1888 1774
1889 1775 return distx, disty, dist, ang
1890 1776 #Calculo de Matrices
1891 1777 # nPairs = len(pairs)
1892 1778 # ang1 = numpy.zeros((nPairs, 2, 1))
1893 1779 # dist1 = numpy.zeros((nPairs, 2, 1))
1894 1780 #
1895 1781 # for j in range(nPairs):
1896 1782 # dist1[j,0,0] = dist[pairs[j][0]]
1897 1783 # dist1[j,1,0] = dist[pairs[j][1]]
1898 1784 # ang1[j,0,0] = ang[pairs[j][0]]
1899 1785 # ang1[j,1,0] = ang[pairs[j][1]]
1900 1786 #
1901 1787 # return distx,disty, dist1,ang1
1902 1788
1903 1789
1904 1790 def __calculateVelVer(self, phase, lagTRange, _lambda):
1905 1791
1906 1792 Ts = lagTRange[1] - lagTRange[0]
1907 1793 velW = -_lambda*phase/(4*math.pi*Ts)
1908 1794
1909 1795 return velW
1910 1796
1911 1797 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1912 1798 nPairs = tau1.shape[0]
1913 1799 nHeights = tau1.shape[1]
1914 1800 vel = numpy.zeros((nPairs,3,nHeights))
1915 1801 dist1 = numpy.reshape(dist, (dist.size,1))
1916 1802
1917 1803 angCos = numpy.cos(ang)
1918 1804 angSin = numpy.sin(ang)
1919 1805
1920 1806 vel0 = dist1*tau1/(2*tau2**2)
1921 1807 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1922 1808 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1923 1809
1924 1810 ind = numpy.where(numpy.isinf(vel))
1925 1811 vel[ind] = numpy.nan
1926 1812
1927 1813 return vel
1928 1814
1929 1815 # def __getPairsAutoCorr(self, pairsList, nChannels):
1930 1816 #
1931 1817 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1932 1818 #
1933 1819 # for l in range(len(pairsList)):
1934 1820 # firstChannel = pairsList[l][0]
1935 1821 # secondChannel = pairsList[l][1]
1936 1822 #
1937 1823 # #Obteniendo pares de Autocorrelacion
1938 1824 # if firstChannel == secondChannel:
1939 1825 # pairsAutoCorr[firstChannel] = int(l)
1940 1826 #
1941 1827 # pairsAutoCorr = pairsAutoCorr.astype(int)
1942 1828 #
1943 1829 # pairsCrossCorr = range(len(pairsList))
1944 1830 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1945 1831 #
1946 1832 # return pairsAutoCorr, pairsCrossCorr
1947 1833
1948 1834 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1949 1835 def techniqueSA(self, kwargs):
1950 1836
1951 1837 """
1952 1838 Function that implements Spaced Antenna (SA) technique.
1953 1839
1954 1840 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1955 1841 Direction correction (if necessary), Ranges and SNR
1956 1842
1957 1843 Output: Winds estimation (Zonal, Meridional and Vertical)
1958 1844
1959 1845 Parameters affected: Winds
1960 1846 """
1961 1847 position_x = kwargs['positionX']
1962 1848 position_y = kwargs['positionY']
1963 1849 azimuth = kwargs['azimuth']
1964 1850
1965 1851 if 'correctFactor' in kwargs:
1966 1852 correctFactor = kwargs['correctFactor']
1967 1853 else:
1968 1854 correctFactor = 1
1969 1855
1970 1856 groupList = kwargs['groupList']
1971 1857 pairs_ccf = groupList[1]
1972 1858 tau = kwargs['tau']
1973 1859 _lambda = kwargs['_lambda']
1974 1860
1975 1861 #Cross Correlation pairs obtained
1976 1862 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
1977 1863 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1978 1864 # pairsSelArray = numpy.array(pairsSelected)
1979 1865 # pairs = []
1980 1866 #
1981 1867 # #Wind estimation pairs obtained
1982 1868 # for i in range(pairsSelArray.shape[0]/2):
1983 1869 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1984 1870 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1985 1871 # pairs.append((ind1,ind2))
1986 1872
1987 1873 indtau = tau.shape[0]/2
1988 1874 tau1 = tau[:indtau,:]
1989 1875 tau2 = tau[indtau:-1,:]
1990 1876 # tau1 = tau1[pairs,:]
1991 1877 # tau2 = tau2[pairs,:]
1992 1878 phase1 = tau[-1,:]
1993 1879
1994 1880 #---------------------------------------------------------------------
1995 1881 #Metodo Directo
1996 1882 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
1997 1883 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1998 1884 winds = stats.nanmean(winds, axis=0)
1999 1885 #---------------------------------------------------------------------
2000 1886 #Metodo General
2001 1887 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
2002 1888 # #Calculo Coeficientes de Funcion de Correlacion
2003 1889 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
2004 1890 # #Calculo de Velocidades
2005 1891 # winds = self.calculateVelUV(F,G,A,B,H)
2006 1892
2007 1893 #---------------------------------------------------------------------
2008 1894 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
2009 1895 winds = correctFactor*winds
2010 1896 return winds
2011 1897
2012 1898 def __checkTime(self, currentTime, paramInterval, outputInterval):
2013 1899
2014 1900 dataTime = currentTime + paramInterval
2015 1901 deltaTime = dataTime - self.__initime
2016 1902
2017 1903 if deltaTime >= outputInterval or deltaTime < 0:
2018 1904 self.__dataReady = True
2019 1905 return
2020 1906
2021 1907 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
2022 1908 '''
2023 1909 Function that implements winds estimation technique with detected meteors.
2024 1910
2025 1911 Input: Detected meteors, Minimum meteor quantity to wind estimation
2026 1912
2027 1913 Output: Winds estimation (Zonal and Meridional)
2028 1914
2029 1915 Parameters affected: Winds
2030 1916 '''
2031 1917 #Settings
2032 1918 nInt = (heightMax - heightMin)/2
2033 1919 nInt = int(nInt)
2034 1920 winds = numpy.zeros((2,nInt))*numpy.nan
2035 1921
2036 1922 #Filter errors
2037 1923 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
2038 1924 finalMeteor = arrayMeteor[error,:]
2039 1925
2040 1926 #Meteor Histogram
2041 1927 finalHeights = finalMeteor[:,2]
2042 1928 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
2043 1929 nMeteorsPerI = hist[0]
2044 1930 heightPerI = hist[1]
2045 1931
2046 1932 #Sort of meteors
2047 1933 indSort = finalHeights.argsort()
2048 1934 finalMeteor2 = finalMeteor[indSort,:]
2049 1935
2050 1936 # Calculating winds
2051 1937 ind1 = 0
2052 1938 ind2 = 0
2053 1939
2054 1940 for i in range(nInt):
2055 1941 nMet = nMeteorsPerI[i]
2056 1942 ind1 = ind2
2057 1943 ind2 = ind1 + nMet
2058 1944
2059 1945 meteorAux = finalMeteor2[ind1:ind2,:]
2060 1946
2061 1947 if meteorAux.shape[0] >= meteorThresh:
2062 1948 vel = meteorAux[:, 6]
2063 1949 zen = meteorAux[:, 4]*numpy.pi/180
2064 1950 azim = meteorAux[:, 3]*numpy.pi/180
2065 1951
2066 1952 n = numpy.cos(zen)
2067 1953 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
2068 1954 # l = m*numpy.tan(azim)
2069 1955 l = numpy.sin(zen)*numpy.sin(azim)
2070 1956 m = numpy.sin(zen)*numpy.cos(azim)
2071 1957
2072 1958 A = numpy.vstack((l, m)).transpose()
2073 1959 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
2074 1960 windsAux = numpy.dot(A1, vel)
2075 1961
2076 1962 winds[0,i] = windsAux[0]
2077 1963 winds[1,i] = windsAux[1]
2078 1964
2079 1965 return winds, heightPerI[:-1]
2080 1966
2081 1967 def techniqueNSM_SA(self, **kwargs):
2082 1968 metArray = kwargs['metArray']
2083 1969 heightList = kwargs['heightList']
2084 1970 timeList = kwargs['timeList']
2085 1971
2086 1972 rx_location = kwargs['rx_location']
2087 1973 groupList = kwargs['groupList']
2088 1974 azimuth = kwargs['azimuth']
2089 1975 dfactor = kwargs['dfactor']
2090 1976 k = kwargs['k']
2091 1977
2092 1978 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
2093 1979 d = dist*dfactor
2094 1980 #Phase calculation
2095 1981 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
2096 1982
2097 1983 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
2098 1984
2099 1985 velEst = numpy.zeros((heightList.size,2))*numpy.nan
2100 1986 azimuth1 = azimuth1*numpy.pi/180
2101 1987
2102 1988 for i in range(heightList.size):
2103 1989 h = heightList[i]
2104 1990 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
2105 1991 metHeight = metArray1[indH,:]
2106 1992 if metHeight.shape[0] >= 2:
2107 1993 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
2108 1994 iazim = metHeight[:,1].astype(int)
2109 1995 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
2110 1996 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
2111 1997 A = numpy.asmatrix(A)
2112 1998 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
2113 1999 velHor = numpy.dot(A1,velAux)
2114 2000
2115 2001 velEst[i,:] = numpy.squeeze(velHor)
2116 2002 return velEst
2117 2003
2118 2004 def __getPhaseSlope(self, metArray, heightList, timeList):
2119 2005 meteorList = []
2120 2006 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
2121 2007 #Putting back together the meteor matrix
2122 2008 utctime = metArray[:,0]
2123 2009 uniqueTime = numpy.unique(utctime)
2124 2010
2125 2011 phaseDerThresh = 0.5
2126 2012 ippSeconds = timeList[1] - timeList[0]
2127 2013 sec = numpy.where(timeList>1)[0][0]
2128 2014 nPairs = metArray.shape[1] - 6
2129 2015 nHeights = len(heightList)
2130 2016
2131 2017 for t in uniqueTime:
2132 2018 metArray1 = metArray[utctime==t,:]
2133 2019 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
2134 2020 tmet = metArray1[:,1].astype(int)
2135 2021 hmet = metArray1[:,2].astype(int)
2136 2022
2137 2023 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
2138 2024 metPhase[:,:] = numpy.nan
2139 2025 metPhase[:,hmet,tmet] = metArray1[:,6:].T
2140 2026
2141 2027 #Delete short trails
2142 2028 metBool = ~numpy.isnan(metPhase[0,:,:])
2143 2029 heightVect = numpy.sum(metBool, axis = 1)
2144 2030 metBool[heightVect<sec,:] = False
2145 2031 metPhase[:,heightVect<sec,:] = numpy.nan
2146 2032
2147 2033 #Derivative
2148 2034 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
2149 2035 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
2150 2036 metPhase[phDerAux] = numpy.nan
2151 2037
2152 2038 #--------------------------METEOR DETECTION -----------------------------------------
2153 2039 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
2154 2040
2155 2041 for p in numpy.arange(nPairs):
2156 2042 phase = metPhase[p,:,:]
2157 2043 phDer = metDer[p,:,:]
2158 2044
2159 2045 for h in indMet:
2160 2046 height = heightList[h]
2161 2047 phase1 = phase[h,:] #82
2162 2048 phDer1 = phDer[h,:]
2163 2049
2164 2050 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
2165 2051
2166 2052 indValid = numpy.where(~numpy.isnan(phase1))[0]
2167 2053 initMet = indValid[0]
2168 2054 endMet = 0
2169 2055
2170 2056 for i in range(len(indValid)-1):
2171 2057
2172 2058 #Time difference
2173 2059 inow = indValid[i]
2174 2060 inext = indValid[i+1]
2175 2061 idiff = inext - inow
2176 2062 #Phase difference
2177 2063 phDiff = numpy.abs(phase1[inext] - phase1[inow])
2178 2064
2179 2065 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
2180 2066 sizeTrail = inow - initMet + 1
2181 2067 if sizeTrail>3*sec: #Too short meteors
2182 2068 x = numpy.arange(initMet,inow+1)*ippSeconds
2183 2069 y = phase1[initMet:inow+1]
2184 2070 ynnan = ~numpy.isnan(y)
2185 2071 x = x[ynnan]
2186 2072 y = y[ynnan]
2187 2073 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
2188 2074 ylin = x*slope + intercept
2189 2075 rsq = r_value**2
2190 2076 if rsq > 0.5:
2191 2077 vel = slope#*height*1000/(k*d)
2192 2078 estAux = numpy.array([utctime,p,height, vel, rsq])
2193 2079 meteorList.append(estAux)
2194 2080 initMet = inext
2195 2081 metArray2 = numpy.array(meteorList)
2196 2082
2197 2083 return metArray2
2198 2084
2199 2085 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
2200 2086
2201 2087 azimuth1 = numpy.zeros(len(pairslist))
2202 2088 dist = numpy.zeros(len(pairslist))
2203 2089
2204 2090 for i in range(len(rx_location)):
2205 2091 ch0 = pairslist[i][0]
2206 2092 ch1 = pairslist[i][1]
2207 2093
2208 2094 diffX = rx_location[ch0][0] - rx_location[ch1][0]
2209 2095 diffY = rx_location[ch0][1] - rx_location[ch1][1]
2210 2096 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
2211 2097 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
2212 2098
2213 2099 azimuth1 -= azimuth0
2214 2100 return azimuth1, dist
2215 2101
2216 2102 def techniqueNSM_DBS(self, **kwargs):
2217 2103 metArray = kwargs['metArray']
2218 2104 heightList = kwargs['heightList']
2219 2105 timeList = kwargs['timeList']
2220 2106 azimuth = kwargs['azimuth']
2221 2107 theta_x = numpy.array(kwargs['theta_x'])
2222 2108 theta_y = numpy.array(kwargs['theta_y'])
2223 2109
2224 2110 utctime = metArray[:,0]
2225 2111 cmet = metArray[:,1].astype(int)
2226 2112 hmet = metArray[:,3].astype(int)
2227 2113 SNRmet = metArray[:,4]
2228 2114 vmet = metArray[:,5]
2229 2115 spcmet = metArray[:,6]
2230 2116
2231 2117 nChan = numpy.max(cmet) + 1
2232 2118 nHeights = len(heightList)
2233 2119
2234 2120 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
2235 2121 hmet = heightList[hmet]
2236 2122 h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
2237 2123
2238 2124 velEst = numpy.zeros((heightList.size,2))*numpy.nan
2239 2125
2240 2126 for i in range(nHeights - 1):
2241 2127 hmin = heightList[i]
2242 2128 hmax = heightList[i + 1]
2243 2129
2244 2130 thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
2245 2131 indthisH = numpy.where(thisH)
2246 2132
2247 2133 if numpy.size(indthisH) > 3:
2248 2134
2249 2135 vel_aux = vmet[thisH]
2250 2136 chan_aux = cmet[thisH]
2251 2137 cosu_aux = dir_cosu[chan_aux]
2252 2138 cosv_aux = dir_cosv[chan_aux]
2253 2139 cosw_aux = dir_cosw[chan_aux]
2254 2140
2255 2141 nch = numpy.size(numpy.unique(chan_aux))
2256 2142 if nch > 1:
2257 2143 A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
2258 2144 velEst[i,:] = numpy.dot(A,vel_aux)
2259 2145
2260 2146 return velEst
2261 2147
2262 2148 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
2263 2149
2264 2150 param = dataOut.data_param
2265 2151 if dataOut.abscissaList != None:
2266 2152 absc = dataOut.abscissaList[:-1]
2267 2153 # noise = dataOut.noise
2268 2154 heightList = dataOut.heightList
2269 2155 SNR = dataOut.data_snr
2270 2156
2271 2157 if technique == 'DBS':
2272 2158
2273 2159 kwargs['velRadial'] = param[:,1,:] #Radial velocity
2274 2160 kwargs['heightList'] = heightList
2275 2161 kwargs['SNR'] = SNR
2276 2162
2277 2163 dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
2278 2164 dataOut.utctimeInit = dataOut.utctime
2279 2165 dataOut.outputInterval = dataOut.paramInterval
2280 2166
2281 2167 elif technique == 'SA':
2282 2168
2283 2169 #Parameters
2284 2170 # position_x = kwargs['positionX']
2285 2171 # position_y = kwargs['positionY']
2286 2172 # azimuth = kwargs['azimuth']
2287 2173 #
2288 2174 # if kwargs.has_key('crosspairsList'):
2289 2175 # pairs = kwargs['crosspairsList']
2290 2176 # else:
2291 2177 # pairs = None
2292 2178 #
2293 2179 # if kwargs.has_key('correctFactor'):
2294 2180 # correctFactor = kwargs['correctFactor']
2295 2181 # else:
2296 2182 # correctFactor = 1
2297 2183
2298 2184 # tau = dataOut.data_param
2299 2185 # _lambda = dataOut.C/dataOut.frequency
2300 2186 # pairsList = dataOut.groupList
2301 2187 # nChannels = dataOut.nChannels
2302 2188
2303 2189 kwargs['groupList'] = dataOut.groupList
2304 2190 kwargs['tau'] = dataOut.data_param
2305 2191 kwargs['_lambda'] = dataOut.C/dataOut.frequency
2306 2192 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
2307 2193 dataOut.data_output = self.techniqueSA(kwargs)
2308 2194 dataOut.utctimeInit = dataOut.utctime
2309 2195 dataOut.outputInterval = dataOut.timeInterval
2310 2196
2311 2197 elif technique == 'Meteors':
2312 2198 dataOut.flagNoData = True
2313 2199 self.__dataReady = False
2314 2200
2315 2201 if 'nHours' in kwargs:
2316 2202 nHours = kwargs['nHours']
2317 2203 else:
2318 2204 nHours = 1
2319 2205
2320 2206 if 'meteorsPerBin' in kwargs:
2321 2207 meteorThresh = kwargs['meteorsPerBin']
2322 2208 else:
2323 2209 meteorThresh = 6
2324 2210
2325 2211 if 'hmin' in kwargs:
2326 2212 hmin = kwargs['hmin']
2327 2213 else: hmin = 70
2328 2214 if 'hmax' in kwargs:
2329 2215 hmax = kwargs['hmax']
2330 2216 else: hmax = 110
2331 2217
2332 2218 dataOut.outputInterval = nHours*3600
2333 2219
2334 2220 if self.__isConfig == False:
2335 2221 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2336 2222 #Get Initial LTC time
2337 2223 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2338 2224 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2339 2225
2340 2226 self.__isConfig = True
2341 2227
2342 2228 if self.__buffer is None:
2343 2229 self.__buffer = dataOut.data_param
2344 2230 self.__firstdata = copy.copy(dataOut)
2345 2231
2346 2232 else:
2347 2233 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2348 2234
2349 2235 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2350 2236
2351 2237 if self.__dataReady:
2352 2238 dataOut.utctimeInit = self.__initime
2353 2239
2354 2240 self.__initime += dataOut.outputInterval #to erase time offset
2355 2241
2356 2242 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
2357 2243 dataOut.flagNoData = False
2358 2244 self.__buffer = None
2359 2245
2360 2246 elif technique == 'Meteors1':
2361 2247 dataOut.flagNoData = True
2362 2248 self.__dataReady = False
2363 2249
2364 2250 if 'nMins' in kwargs:
2365 2251 nMins = kwargs['nMins']
2366 2252 else: nMins = 20
2367 2253 if 'rx_location' in kwargs:
2368 2254 rx_location = kwargs['rx_location']
2369 2255 else: rx_location = [(0,1),(1,1),(1,0)]
2370 2256 if 'azimuth' in kwargs:
2371 2257 azimuth = kwargs['azimuth']
2372 2258 else: azimuth = 51.06
2373 2259 if 'dfactor' in kwargs:
2374 2260 dfactor = kwargs['dfactor']
2375 2261 if 'mode' in kwargs:
2376 2262 mode = kwargs['mode']
2377 2263 if 'theta_x' in kwargs:
2378 2264 theta_x = kwargs['theta_x']
2379 2265 if 'theta_y' in kwargs:
2380 2266 theta_y = kwargs['theta_y']
2381 2267 else: mode = 'SA'
2382 2268
2383 2269 #Borrar luego esto
2384 2270 if dataOut.groupList is None:
2385 2271 dataOut.groupList = [(0,1),(0,2),(1,2)]
2386 2272 groupList = dataOut.groupList
2387 2273 C = 3e8
2388 2274 freq = 50e6
2389 2275 lamb = C/freq
2390 2276 k = 2*numpy.pi/lamb
2391 2277
2392 2278 timeList = dataOut.abscissaList
2393 2279 heightList = dataOut.heightList
2394 2280
2395 2281 if self.__isConfig == False:
2396 2282 dataOut.outputInterval = nMins*60
2397 2283 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2398 2284 #Get Initial LTC time
2399 2285 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2400 2286 minuteAux = initime.minute
2401 2287 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
2402 2288 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2403 2289
2404 2290 self.__isConfig = True
2405 2291
2406 2292 if self.__buffer is None:
2407 2293 self.__buffer = dataOut.data_param
2408 2294 self.__firstdata = copy.copy(dataOut)
2409 2295
2410 2296 else:
2411 2297 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2412 2298
2413 2299 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2414 2300
2415 2301 if self.__dataReady:
2416 2302 dataOut.utctimeInit = self.__initime
2417 2303 self.__initime += dataOut.outputInterval #to erase time offset
2418 2304
2419 2305 metArray = self.__buffer
2420 2306 if mode == 'SA':
2421 2307 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
2422 2308 elif mode == 'DBS':
2423 2309 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
2424 2310 dataOut.data_output = dataOut.data_output.T
2425 2311 dataOut.flagNoData = False
2426 2312 self.__buffer = None
2427 2313
2428 2314 return
2429 2315
2430 2316 class EWDriftsEstimation(Operation):
2431 2317
2432 2318 def __init__(self):
2433 2319 Operation.__init__(self)
2434 2320
2435 2321 def __correctValues(self, heiRang, phi, velRadial, SNR):
2436 2322 listPhi = phi.tolist()
2437 2323 maxid = listPhi.index(max(listPhi))
2438 2324 minid = listPhi.index(min(listPhi))
2439 2325
2440 2326 rango = list(range(len(phi)))
2441 2327 # rango = numpy.delete(rango,maxid)
2442 2328
2443 2329 heiRang1 = heiRang*math.cos(phi[maxid])
2444 2330 heiRangAux = heiRang*math.cos(phi[minid])
2445 2331 indOut = (heiRang1 < heiRangAux[0]).nonzero()
2446 2332 heiRang1 = numpy.delete(heiRang1,indOut)
2447 2333
2448 2334 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2449 2335 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2450 2336
2451 2337 for i in rango:
2452 2338 x = heiRang*math.cos(phi[i])
2453 2339 y1 = velRadial[i,:]
2454 2340 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2455 2341
2456 2342 x1 = heiRang1
2457 2343 y11 = f1(x1)
2458 2344
2459 2345 y2 = SNR[i,:]
2460 2346 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2461 2347 y21 = f2(x1)
2462 2348
2463 2349 velRadial1[i,:] = y11
2464 2350 SNR1[i,:] = y21
2465 2351
2466 2352 return heiRang1, velRadial1, SNR1
2467 2353
2468 2354 def run(self, dataOut, zenith, zenithCorrection):
2469 2355 heiRang = dataOut.heightList
2470 2356 velRadial = dataOut.data_param[:,3,:]
2471 2357 SNR = dataOut.data_snr
2472 2358
2473 2359 zenith = numpy.array(zenith)
2474 2360 zenith -= zenithCorrection
2475 2361 zenith *= numpy.pi/180
2476 2362
2477 2363 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2478 2364
2479 2365 alp = zenith[0]
2480 2366 bet = zenith[1]
2481 2367
2482 2368 w_w = velRadial1[0,:]
2483 2369 w_e = velRadial1[1,:]
2484 2370
2485 2371 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
2486 2372 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
2487 2373
2488 2374 winds = numpy.vstack((u,w))
2489 2375
2490 2376 dataOut.heightList = heiRang1
2491 2377 dataOut.data_output = winds
2492 2378 dataOut.data_snr = SNR1
2493 2379
2494 2380 dataOut.utctimeInit = dataOut.utctime
2495 2381 dataOut.outputInterval = dataOut.timeInterval
2496 2382 return
2497 2383
2498 2384 #--------------- Non Specular Meteor ----------------
2499 2385
2500 2386 class NonSpecularMeteorDetection(Operation):
2501 2387
2502 2388 def run(self, dataOut, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
2503 2389 data_acf = dataOut.data_pre[0]
2504 2390 data_ccf = dataOut.data_pre[1]
2505 2391 pairsList = dataOut.groupList[1]
2506 2392
2507 2393 lamb = dataOut.C/dataOut.frequency
2508 2394 tSamp = dataOut.ippSeconds*dataOut.nCohInt
2509 2395 paramInterval = dataOut.paramInterval
2510 2396
2511 2397 nChannels = data_acf.shape[0]
2512 2398 nLags = data_acf.shape[1]
2513 2399 nProfiles = data_acf.shape[2]
2514 2400 nHeights = dataOut.nHeights
2515 2401 nCohInt = dataOut.nCohInt
2516 2402 sec = numpy.round(nProfiles/dataOut.paramInterval)
2517 2403 heightList = dataOut.heightList
2518 2404 ippSeconds = dataOut.ippSeconds*dataOut.nCohInt*dataOut.nAvg
2519 2405 utctime = dataOut.utctime
2520 2406
2521 2407 dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
2522 2408
2523 2409 #------------------------ SNR --------------------------------------
2524 2410 power = data_acf[:,0,:,:].real
2525 2411 noise = numpy.zeros(nChannels)
2526 2412 SNR = numpy.zeros(power.shape)
2527 2413 for i in range(nChannels):
2528 2414 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
2529 2415 SNR[i] = (power[i]-noise[i])/noise[i]
2530 2416 SNRm = numpy.nanmean(SNR, axis = 0)
2531 2417 SNRdB = 10*numpy.log10(SNR)
2532 2418
2533 2419 if mode == 'SA':
2534 2420 dataOut.groupList = dataOut.groupList[1]
2535 2421 nPairs = data_ccf.shape[0]
2536 2422 #---------------------- Coherence and Phase --------------------------
2537 2423 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
2538 2424 # phase1 = numpy.copy(phase)
2539 2425 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
2540 2426
2541 2427 for p in range(nPairs):
2542 2428 ch0 = pairsList[p][0]
2543 2429 ch1 = pairsList[p][1]
2544 2430 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
2545 2431 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
2546 2432 # phase1[p,:,:] = numpy.angle(ccf) #median filter
2547 2433 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
2548 2434 # coh1[p,:,:] = numpy.abs(ccf) #median filter
2549 2435 coh = numpy.nanmax(coh1, axis = 0)
2550 2436 # struc = numpy.ones((5,1))
2551 2437 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
2552 2438 #---------------------- Radial Velocity ----------------------------
2553 2439 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
2554 2440 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
2555 2441
2556 2442 if allData:
2557 2443 boolMetFin = ~numpy.isnan(SNRm)
2558 2444 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2559 2445 else:
2560 2446 #------------------------ Meteor mask ---------------------------------
2561 2447 # #SNR mask
2562 2448 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
2563 2449 #
2564 2450 # #Erase small objects
2565 2451 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
2566 2452 #
2567 2453 # auxEEJ = numpy.sum(boolMet1,axis=0)
2568 2454 # indOver = auxEEJ>nProfiles*0.8 #Use this later
2569 2455 # indEEJ = numpy.where(indOver)[0]
2570 2456 # indNEEJ = numpy.where(~indOver)[0]
2571 2457 #
2572 2458 # boolMetFin = boolMet1
2573 2459 #
2574 2460 # if indEEJ.size > 0:
2575 2461 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
2576 2462 #
2577 2463 # boolMet2 = coh > cohThresh
2578 2464 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
2579 2465 #
2580 2466 # #Final Meteor mask
2581 2467 # boolMetFin = boolMet1|boolMet2
2582 2468
2583 2469 #Coherence mask
2584 2470 boolMet1 = coh > 0.75
2585 2471 struc = numpy.ones((30,1))
2586 2472 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
2587 2473
2588 2474 #Derivative mask
2589 2475 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2590 2476 boolMet2 = derPhase < 0.2
2591 2477 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
2592 2478 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
2593 2479 boolMet2 = ndimage.median_filter(boolMet2,size=5)
2594 2480 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
2595 2481 # #Final mask
2596 2482 # boolMetFin = boolMet2
2597 2483 boolMetFin = boolMet1&boolMet2
2598 2484 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
2599 2485 #Creating data_param
2600 2486 coordMet = numpy.where(boolMetFin)
2601 2487
2602 2488 tmet = coordMet[0]
2603 2489 hmet = coordMet[1]
2604 2490
2605 2491 data_param = numpy.zeros((tmet.size, 6 + nPairs))
2606 2492 data_param[:,0] = utctime
2607 2493 data_param[:,1] = tmet
2608 2494 data_param[:,2] = hmet
2609 2495 data_param[:,3] = SNRm[tmet,hmet]
2610 2496 data_param[:,4] = velRad[tmet,hmet]
2611 2497 data_param[:,5] = coh[tmet,hmet]
2612 2498 data_param[:,6:] = phase[:,tmet,hmet].T
2613 2499
2614 2500 elif mode == 'DBS':
2615 2501 dataOut.groupList = numpy.arange(nChannels)
2616 2502
2617 2503 #Radial Velocities
2618 2504 phase = numpy.angle(data_acf[:,1,:,:])
2619 2505 # phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
2620 2506 velRad = phase*lamb/(4*numpy.pi*tSamp)
2621 2507
2622 2508 #Spectral width
2623 2509 # acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
2624 2510 # acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
2625 2511 acf1 = data_acf[:,1,:,:]
2626 2512 acf2 = data_acf[:,2,:,:]
2627 2513
2628 2514 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
2629 2515 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
2630 2516 if allData:
2631 2517 boolMetFin = ~numpy.isnan(SNRdB)
2632 2518 else:
2633 2519 #SNR
2634 2520 boolMet1 = (SNRdB>SNRthresh) #SNR mask
2635 2521 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
2636 2522
2637 2523 #Radial velocity
2638 2524 boolMet2 = numpy.abs(velRad) < 20
2639 2525 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
2640 2526
2641 2527 #Spectral Width
2642 2528 boolMet3 = spcWidth < 30
2643 2529 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
2644 2530 # boolMetFin = self.__erase_small(boolMet1, 10,5)
2645 2531 boolMetFin = boolMet1&boolMet2&boolMet3
2646 2532
2647 2533 #Creating data_param
2648 2534 coordMet = numpy.where(boolMetFin)
2649 2535
2650 2536 cmet = coordMet[0]
2651 2537 tmet = coordMet[1]
2652 2538 hmet = coordMet[2]
2653 2539
2654 2540 data_param = numpy.zeros((tmet.size, 7))
2655 2541 data_param[:,0] = utctime
2656 2542 data_param[:,1] = cmet
2657 2543 data_param[:,2] = tmet
2658 2544 data_param[:,3] = hmet
2659 2545 data_param[:,4] = SNR[cmet,tmet,hmet].T
2660 2546 data_param[:,5] = velRad[cmet,tmet,hmet].T
2661 2547 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
2662 2548
2663 2549 # self.dataOut.data_param = data_int
2664 2550 if len(data_param) == 0:
2665 2551 dataOut.flagNoData = True
2666 2552 else:
2667 2553 dataOut.data_param = data_param
2668 2554
2669 2555 def __erase_small(self, binArray, threshX, threshY):
2670 2556 labarray, numfeat = ndimage.measurements.label(binArray)
2671 2557 binArray1 = numpy.copy(binArray)
2672 2558
2673 2559 for i in range(1,numfeat + 1):
2674 2560 auxBin = (labarray==i)
2675 2561 auxSize = auxBin.sum()
2676 2562
2677 2563 x,y = numpy.where(auxBin)
2678 2564 widthX = x.max() - x.min()
2679 2565 widthY = y.max() - y.min()
2680 2566
2681 2567 #width X: 3 seg -> 12.5*3
2682 2568 #width Y:
2683 2569
2684 2570 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
2685 2571 binArray1[auxBin] = False
2686 2572
2687 2573 return binArray1
2688 2574
2689 2575 #--------------- Specular Meteor ----------------
2690 2576
2691 2577 class SMDetection(Operation):
2692 2578 '''
2693 2579 Function DetectMeteors()
2694 2580 Project developed with paper:
2695 2581 HOLDSWORTH ET AL. 2004
2696 2582
2697 2583 Input:
2698 2584 self.dataOut.data_pre
2699 2585
2700 2586 centerReceiverIndex: From the channels, which is the center receiver
2701 2587
2702 2588 hei_ref: Height reference for the Beacon signal extraction
2703 2589 tauindex:
2704 2590 predefinedPhaseShifts: Predefined phase offset for the voltge signals
2705 2591
2706 2592 cohDetection: Whether to user Coherent detection or not
2707 2593 cohDet_timeStep: Coherent Detection calculation time step
2708 2594 cohDet_thresh: Coherent Detection phase threshold to correct phases
2709 2595
2710 2596 noise_timeStep: Noise calculation time step
2711 2597 noise_multiple: Noise multiple to define signal threshold
2712 2598
2713 2599 multDet_timeLimit: Multiple Detection Removal time limit in seconds
2714 2600 multDet_rangeLimit: Multiple Detection Removal range limit in km
2715 2601
2716 2602 phaseThresh: Maximum phase difference between receiver to be consider a meteor
2717 2603 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
2718 2604
2719 2605 hmin: Minimum Height of the meteor to use it in the further wind estimations
2720 2606 hmax: Maximum Height of the meteor to use it in the further wind estimations
2721 2607 azimuth: Azimuth angle correction
2722 2608
2723 2609 Affected:
2724 2610 self.dataOut.data_param
2725 2611
2726 2612 Rejection Criteria (Errors):
2727 2613 0: No error; analysis OK
2728 2614 1: SNR < SNR threshold
2729 2615 2: angle of arrival (AOA) ambiguously determined
2730 2616 3: AOA estimate not feasible
2731 2617 4: Large difference in AOAs obtained from different antenna baselines
2732 2618 5: echo at start or end of time series
2733 2619 6: echo less than 5 examples long; too short for analysis
2734 2620 7: echo rise exceeds 0.3s
2735 2621 8: echo decay time less than twice rise time
2736 2622 9: large power level before echo
2737 2623 10: large power level after echo
2738 2624 11: poor fit to amplitude for estimation of decay time
2739 2625 12: poor fit to CCF phase variation for estimation of radial drift velocity
2740 2626 13: height unresolvable echo: not valid height within 70 to 110 km
2741 2627 14: height ambiguous echo: more then one possible height within 70 to 110 km
2742 2628 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
2743 2629 16: oscilatory echo, indicating event most likely not an underdense echo
2744 2630
2745 2631 17: phase difference in meteor Reestimation
2746 2632
2747 2633 Data Storage:
2748 2634 Meteors for Wind Estimation (8):
2749 2635 Utc Time | Range Height
2750 2636 Azimuth Zenith errorCosDir
2751 2637 VelRad errorVelRad
2752 2638 Phase0 Phase1 Phase2 Phase3
2753 2639 TypeError
2754 2640
2755 2641 '''
2756 2642
2757 2643 def run(self, dataOut, hei_ref = None, tauindex = 0,
2758 2644 phaseOffsets = None,
2759 2645 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
2760 2646 noise_timeStep = 4, noise_multiple = 4,
2761 2647 multDet_timeLimit = 1, multDet_rangeLimit = 3,
2762 2648 phaseThresh = 20, SNRThresh = 5,
2763 2649 hmin = 50, hmax=150, azimuth = 0,
2764 2650 channelPositions = None) :
2765 2651
2766 2652
2767 2653 #Getting Pairslist
2768 2654 if channelPositions is None:
2769 2655 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2770 2656 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2771 2657 meteorOps = SMOperations()
2772 2658 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2773 2659 heiRang = dataOut.heightList
2774 2660 #Get Beacon signal - No Beacon signal anymore
2775 2661 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
2776 2662 #
2777 2663 # if hei_ref != None:
2778 2664 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
2779 2665 #
2780 2666
2781 2667
2782 2668 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
2783 2669 # see if the user put in pre defined phase shifts
2784 2670 voltsPShift = dataOut.data_pre.copy()
2785 2671
2786 2672 # if predefinedPhaseShifts != None:
2787 2673 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
2788 2674 #
2789 2675 # # elif beaconPhaseShifts:
2790 2676 # # #get hardware phase shifts using beacon signal
2791 2677 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
2792 2678 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
2793 2679 #
2794 2680 # else:
2795 2681 # hardwarePhaseShifts = numpy.zeros(5)
2796 2682 #
2797 2683 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
2798 2684 # for i in range(self.dataOut.data_pre.shape[0]):
2799 2685 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
2800 2686
2801 2687 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
2802 2688
2803 2689 #Remove DC
2804 2690 voltsDC = numpy.mean(voltsPShift,1)
2805 2691 voltsDC = numpy.mean(voltsDC,1)
2806 2692 for i in range(voltsDC.shape[0]):
2807 2693 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
2808 2694
2809 2695 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
2810 2696 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
2811 2697
2812 2698 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
2813 2699 #Coherent Detection
2814 2700 if cohDetection:
2815 2701 #use coherent detection to get the net power
2816 2702 cohDet_thresh = cohDet_thresh*numpy.pi/180
2817 2703 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
2818 2704
2819 2705 #Non-coherent detection!
2820 2706 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
2821 2707 #********** END OF COH/NON-COH POWER CALCULATION**********************
2822 2708
2823 2709 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
2824 2710 #Get noise
2825 2711 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
2826 2712 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
2827 2713 #Get signal threshold
2828 2714 signalThresh = noise_multiple*noise
2829 2715 #Meteor echoes detection
2830 2716 listMeteors = self.__findMeteors(powerNet, signalThresh)
2831 2717 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
2832 2718
2833 2719 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
2834 2720 #Parameters
2835 2721 heiRange = dataOut.heightList
2836 2722 rangeInterval = heiRange[1] - heiRange[0]
2837 2723 rangeLimit = multDet_rangeLimit/rangeInterval
2838 2724 timeLimit = multDet_timeLimit/dataOut.timeInterval
2839 2725 #Multiple detection removals
2840 2726 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
2841 2727 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
2842 2728
2843 2729 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
2844 2730 #Parameters
2845 2731 phaseThresh = phaseThresh*numpy.pi/180
2846 2732 thresh = [phaseThresh, noise_multiple, SNRThresh]
2847 2733 #Meteor reestimation (Errors N 1, 6, 12, 17)
2848 2734 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
2849 2735 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
2850 2736 #Estimation of decay times (Errors N 7, 8, 11)
2851 2737 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
2852 2738 #******************* END OF METEOR REESTIMATION *******************
2853 2739
2854 2740 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
2855 2741 #Calculating Radial Velocity (Error N 15)
2856 2742 radialStdThresh = 10
2857 2743 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
2858 2744
2859 2745 if len(listMeteors4) > 0:
2860 2746 #Setting New Array
2861 2747 date = dataOut.utctime
2862 2748 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
2863 2749
2864 2750 #Correcting phase offset
2865 2751 if phaseOffsets != None:
2866 2752 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2867 2753 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2868 2754
2869 2755 #Second Pairslist
2870 2756 pairsList = []
2871 2757 pairx = (0,1)
2872 2758 pairy = (2,3)
2873 2759 pairsList.append(pairx)
2874 2760 pairsList.append(pairy)
2875 2761
2876 2762 jph = numpy.array([0,0,0,0])
2877 2763 h = (hmin,hmax)
2878 2764 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2879 2765
2880 2766 # #Calculate AOA (Error N 3, 4)
2881 2767 # #JONES ET AL. 1998
2882 2768 # error = arrayParameters[:,-1]
2883 2769 # AOAthresh = numpy.pi/8
2884 2770 # phases = -arrayParameters[:,9:13]
2885 2771 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
2886 2772 #
2887 2773 # #Calculate Heights (Error N 13 and 14)
2888 2774 # error = arrayParameters[:,-1]
2889 2775 # Ranges = arrayParameters[:,2]
2890 2776 # zenith = arrayParameters[:,5]
2891 2777 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
2892 2778 # error = arrayParameters[:,-1]
2893 2779 #********************* END OF PARAMETERS CALCULATION **************************
2894 2780
2895 2781 #***************************+ PASS DATA TO NEXT STEP **********************
2896 2782 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
2897 2783 dataOut.data_param = arrayParameters
2898 2784
2899 2785 if arrayParameters is None:
2900 2786 dataOut.flagNoData = True
2901 2787 else:
2902 2788 dataOut.flagNoData = True
2903 2789
2904 2790 return
2905 2791
2906 2792 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
2907 2793
2908 2794 minIndex = min(newheis[0])
2909 2795 maxIndex = max(newheis[0])
2910 2796
2911 2797 voltage = voltage0[:,:,minIndex:maxIndex+1]
2912 2798 nLength = voltage.shape[1]/n
2913 2799 nMin = 0
2914 2800 nMax = 0
2915 2801 phaseOffset = numpy.zeros((len(pairslist),n))
2916 2802
2917 2803 for i in range(n):
2918 2804 nMax += nLength
2919 2805 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
2920 2806 phaseCCF = numpy.mean(phaseCCF, axis = 2)
2921 2807 phaseOffset[:,i] = phaseCCF.transpose()
2922 2808 nMin = nMax
2923 2809 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
2924 2810
2925 2811 #Remove Outliers
2926 2812 factor = 2
2927 2813 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
2928 2814 dw = numpy.std(wt,axis = 1)
2929 2815 dw = dw.reshape((dw.size,1))
2930 2816 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
2931 2817 phaseOffset[ind] = numpy.nan
2932 2818 phaseOffset = stats.nanmean(phaseOffset, axis=1)
2933 2819
2934 2820 return phaseOffset
2935 2821
2936 2822 def __shiftPhase(self, data, phaseShift):
2937 2823 #this will shift the phase of a complex number
2938 2824 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
2939 2825 return dataShifted
2940 2826
2941 2827 def __estimatePhaseDifference(self, array, pairslist):
2942 2828 nChannel = array.shape[0]
2943 2829 nHeights = array.shape[2]
2944 2830 numPairs = len(pairslist)
2945 2831 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
2946 2832 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
2947 2833
2948 2834 #Correct phases
2949 2835 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
2950 2836 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2951 2837
2952 2838 if indDer[0].shape[0] > 0:
2953 2839 for i in range(indDer[0].shape[0]):
2954 2840 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
2955 2841 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
2956 2842
2957 2843 # for j in range(numSides):
2958 2844 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
2959 2845 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
2960 2846 #
2961 2847 #Linear
2962 2848 phaseInt = numpy.zeros((numPairs,1))
2963 2849 angAllCCF = phaseCCF[:,[0,1,3,4],0]
2964 2850 for j in range(numPairs):
2965 2851 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
2966 2852 phaseInt[j] = fit[1]
2967 2853 #Phase Differences
2968 2854 phaseDiff = phaseInt - phaseCCF[:,2,:]
2969 2855 phaseArrival = phaseInt.reshape(phaseInt.size)
2970 2856
2971 2857 #Dealias
2972 2858 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
2973 2859 # indAlias = numpy.where(phaseArrival > numpy.pi)
2974 2860 # phaseArrival[indAlias] -= 2*numpy.pi
2975 2861 # indAlias = numpy.where(phaseArrival < -numpy.pi)
2976 2862 # phaseArrival[indAlias] += 2*numpy.pi
2977 2863
2978 2864 return phaseDiff, phaseArrival
2979 2865
2980 2866 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
2981 2867 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
2982 2868 #find the phase shifts of each channel over 1 second intervals
2983 2869 #only look at ranges below the beacon signal
2984 2870 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
2985 2871 numBlocks = int(volts.shape[1]/numProfPerBlock)
2986 2872 numHeights = volts.shape[2]
2987 2873 nChannel = volts.shape[0]
2988 2874 voltsCohDet = volts.copy()
2989 2875
2990 2876 pairsarray = numpy.array(pairslist)
2991 2877 indSides = pairsarray[:,1]
2992 2878 # indSides = numpy.array(range(nChannel))
2993 2879 # indSides = numpy.delete(indSides, indCenter)
2994 2880 #
2995 2881 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
2996 2882 listBlocks = numpy.array_split(volts, numBlocks, 1)
2997 2883
2998 2884 startInd = 0
2999 2885 endInd = 0
3000 2886
3001 2887 for i in range(numBlocks):
3002 2888 startInd = endInd
3003 2889 endInd = endInd + listBlocks[i].shape[1]
3004 2890
3005 2891 arrayBlock = listBlocks[i]
3006 2892 # arrayBlockCenter = listCenter[i]
3007 2893
3008 2894 #Estimate the Phase Difference
3009 2895 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
3010 2896 #Phase Difference RMS
3011 2897 arrayPhaseRMS = numpy.abs(phaseDiff)
3012 2898 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
3013 2899 indPhase = numpy.where(phaseRMSaux==4)
3014 2900 #Shifting
3015 2901 if indPhase[0].shape[0] > 0:
3016 2902 for j in range(indSides.size):
3017 2903 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
3018 2904 voltsCohDet[:,startInd:endInd,:] = arrayBlock
3019 2905
3020 2906 return voltsCohDet
3021 2907
3022 2908 def __calculateCCF(self, volts, pairslist ,laglist):
3023 2909
3024 2910 nHeights = volts.shape[2]
3025 2911 nPoints = volts.shape[1]
3026 2912 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
3027 2913
3028 2914 for i in range(len(pairslist)):
3029 2915 volts1 = volts[pairslist[i][0]]
3030 2916 volts2 = volts[pairslist[i][1]]
3031 2917
3032 2918 for t in range(len(laglist)):
3033 2919 idxT = laglist[t]
3034 2920 if idxT >= 0:
3035 2921 vStacked = numpy.vstack((volts2[idxT:,:],
3036 2922 numpy.zeros((idxT, nHeights),dtype='complex')))
3037 2923 else:
3038 2924 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
3039 2925 volts2[:(nPoints + idxT),:]))
3040 2926 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
3041 2927
3042 2928 vStacked = None
3043 2929 return voltsCCF
3044 2930
3045 2931 def __getNoise(self, power, timeSegment, timeInterval):
3046 2932 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
3047 2933 numBlocks = int(power.shape[0]/numProfPerBlock)
3048 2934 numHeights = power.shape[1]
3049 2935
3050 2936 listPower = numpy.array_split(power, numBlocks, 0)
3051 2937 noise = numpy.zeros((power.shape[0], power.shape[1]))
3052 2938 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
3053 2939
3054 2940 startInd = 0
3055 2941 endInd = 0
3056 2942
3057 2943 for i in range(numBlocks): #split por canal
3058 2944 startInd = endInd
3059 2945 endInd = endInd + listPower[i].shape[0]
3060 2946
3061 2947 arrayBlock = listPower[i]
3062 2948 noiseAux = numpy.mean(arrayBlock, 0)
3063 2949 # noiseAux = numpy.median(noiseAux)
3064 2950 # noiseAux = numpy.mean(arrayBlock)
3065 2951 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
3066 2952
3067 2953 noiseAux1 = numpy.mean(arrayBlock)
3068 2954 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
3069 2955
3070 2956 return noise, noise1
3071 2957
3072 2958 def __findMeteors(self, power, thresh):
3073 2959 nProf = power.shape[0]
3074 2960 nHeights = power.shape[1]
3075 2961 listMeteors = []
3076 2962
3077 2963 for i in range(nHeights):
3078 2964 powerAux = power[:,i]
3079 2965 threshAux = thresh[:,i]
3080 2966
3081 2967 indUPthresh = numpy.where(powerAux > threshAux)[0]
3082 2968 indDNthresh = numpy.where(powerAux <= threshAux)[0]
3083 2969
3084 2970 j = 0
3085 2971
3086 2972 while (j < indUPthresh.size - 2):
3087 2973 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
3088 2974 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
3089 2975 indDNthresh = indDNthresh[indDNAux]
3090 2976
3091 2977 if (indDNthresh.size > 0):
3092 2978 indEnd = indDNthresh[0] - 1
3093 2979 indInit = indUPthresh[j]
3094 2980
3095 2981 meteor = powerAux[indInit:indEnd + 1]
3096 2982 indPeak = meteor.argmax() + indInit
3097 2983 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
3098 2984
3099 2985 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
3100 2986 j = numpy.where(indUPthresh == indEnd)[0] + 1
3101 2987 else: j+=1
3102 2988 else: j+=1
3103 2989
3104 2990 return listMeteors
3105 2991
3106 2992 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
3107 2993
3108 2994 arrayMeteors = numpy.asarray(listMeteors)
3109 2995 listMeteors1 = []
3110 2996
3111 2997 while arrayMeteors.shape[0] > 0:
3112 2998 FLAs = arrayMeteors[:,4]
3113 2999 maxFLA = FLAs.argmax()
3114 3000 listMeteors1.append(arrayMeteors[maxFLA,:])
3115 3001
3116 3002 MeteorInitTime = arrayMeteors[maxFLA,1]
3117 3003 MeteorEndTime = arrayMeteors[maxFLA,3]
3118 3004 MeteorHeight = arrayMeteors[maxFLA,0]
3119 3005
3120 3006 #Check neighborhood
3121 3007 maxHeightIndex = MeteorHeight + rangeLimit
3122 3008 minHeightIndex = MeteorHeight - rangeLimit
3123 3009 minTimeIndex = MeteorInitTime - timeLimit
3124 3010 maxTimeIndex = MeteorEndTime + timeLimit
3125 3011
3126 3012 #Check Heights
3127 3013 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
3128 3014 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
3129 3015 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
3130 3016
3131 3017 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
3132 3018
3133 3019 return listMeteors1
3134 3020
3135 3021 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
3136 3022 numHeights = volts.shape[2]
3137 3023 nChannel = volts.shape[0]
3138 3024
3139 3025 thresholdPhase = thresh[0]
3140 3026 thresholdNoise = thresh[1]
3141 3027 thresholdDB = float(thresh[2])
3142 3028
3143 3029 thresholdDB1 = 10**(thresholdDB/10)
3144 3030 pairsarray = numpy.array(pairslist)
3145 3031 indSides = pairsarray[:,1]
3146 3032
3147 3033 pairslist1 = list(pairslist)
3148 3034 pairslist1.append((0,1))
3149 3035 pairslist1.append((3,4))
3150 3036
3151 3037 listMeteors1 = []
3152 3038 listPowerSeries = []
3153 3039 listVoltageSeries = []
3154 3040 #volts has the war data
3155 3041
3156 3042 if frequency == 30e6:
3157 3043 timeLag = 45*10**-3
3158 3044 else:
3159 3045 timeLag = 15*10**-3
3160 3046 lag = numpy.ceil(timeLag/timeInterval)
3161 3047
3162 3048 for i in range(len(listMeteors)):
3163 3049
3164 3050 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
3165 3051 meteorAux = numpy.zeros(16)
3166 3052
3167 3053 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
3168 3054 mHeight = listMeteors[i][0]
3169 3055 mStart = listMeteors[i][1]
3170 3056 mPeak = listMeteors[i][2]
3171 3057 mEnd = listMeteors[i][3]
3172 3058
3173 3059 #get the volt data between the start and end times of the meteor
3174 3060 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
3175 3061 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3176 3062
3177 3063 #3.6. Phase Difference estimation
3178 3064 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
3179 3065
3180 3066 #3.7. Phase difference removal & meteor start, peak and end times reestimated
3181 3067 #meteorVolts0.- all Channels, all Profiles
3182 3068 meteorVolts0 = volts[:,:,mHeight]
3183 3069 meteorThresh = noise[:,mHeight]*thresholdNoise
3184 3070 meteorNoise = noise[:,mHeight]
3185 3071 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
3186 3072 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
3187 3073
3188 3074 #Times reestimation
3189 3075 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
3190 3076 if mStart1.size > 0:
3191 3077 mStart1 = mStart1[-1] + 1
3192 3078
3193 3079 else:
3194 3080 mStart1 = mPeak
3195 3081
3196 3082 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
3197 3083 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
3198 3084 if mEndDecayTime1.size == 0:
3199 3085 mEndDecayTime1 = powerNet0.size
3200 3086 else:
3201 3087 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
3202 3088 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
3203 3089
3204 3090 #meteorVolts1.- all Channels, from start to end
3205 3091 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
3206 3092 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
3207 3093 if meteorVolts2.shape[1] == 0:
3208 3094 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
3209 3095 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
3210 3096 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
3211 3097 ##################### END PARAMETERS REESTIMATION #########################
3212 3098
3213 3099 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
3214 3100 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
3215 3101 if meteorVolts2.shape[1] > 0:
3216 3102 #Phase Difference re-estimation
3217 3103 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
3218 3104 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
3219 3105 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
3220 3106 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
3221 3107 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
3222 3108
3223 3109 #Phase Difference RMS
3224 3110 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
3225 3111 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
3226 3112 #Data from Meteor
3227 3113 mPeak1 = powerNet1.argmax() + mStart1
3228 3114 mPeakPower1 = powerNet1.max()
3229 3115 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
3230 3116 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
3231 3117 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
3232 3118 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
3233 3119 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
3234 3120 #Vectorize
3235 3121 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
3236 3122 meteorAux[7:11] = phaseDiffint[0:4]
3237 3123
3238 3124 #Rejection Criterions
3239 3125 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
3240 3126 meteorAux[-1] = 17
3241 3127 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
3242 3128 meteorAux[-1] = 1
3243 3129
3244 3130
3245 3131 else:
3246 3132 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
3247 3133 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
3248 3134 PowerSeries = 0
3249 3135
3250 3136 listMeteors1.append(meteorAux)
3251 3137 listPowerSeries.append(PowerSeries)
3252 3138 listVoltageSeries.append(meteorVolts1)
3253 3139
3254 3140 return listMeteors1, listPowerSeries, listVoltageSeries
3255 3141
3256 3142 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
3257 3143
3258 3144 threshError = 10
3259 3145 #Depending if it is 30 or 50 MHz
3260 3146 if frequency == 30e6:
3261 3147 timeLag = 45*10**-3
3262 3148 else:
3263 3149 timeLag = 15*10**-3
3264 3150 lag = numpy.ceil(timeLag/timeInterval)
3265 3151
3266 3152 listMeteors1 = []
3267 3153
3268 3154 for i in range(len(listMeteors)):
3269 3155 meteorPower = listPower[i]
3270 3156 meteorAux = listMeteors[i]
3271 3157
3272 3158 if meteorAux[-1] == 0:
3273 3159
3274 3160 try:
3275 3161 indmax = meteorPower.argmax()
3276 3162 indlag = indmax + lag
3277 3163
3278 3164 y = meteorPower[indlag:]
3279 3165 x = numpy.arange(0, y.size)*timeLag
3280 3166
3281 3167 #first guess
3282 3168 a = y[0]
3283 3169 tau = timeLag
3284 3170 #exponential fit
3285 3171 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
3286 3172 y1 = self.__exponential_function(x, *popt)
3287 3173 #error estimation
3288 3174 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
3289 3175
3290 3176 decayTime = popt[1]
3291 3177 riseTime = indmax*timeInterval
3292 3178 meteorAux[11:13] = [decayTime, error]
3293 3179
3294 3180 #Table items 7, 8 and 11
3295 3181 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
3296 3182 meteorAux[-1] = 7
3297 3183 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
3298 3184 meteorAux[-1] = 8
3299 3185 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
3300 3186 meteorAux[-1] = 11
3301 3187
3302 3188
3303 3189 except:
3304 3190 meteorAux[-1] = 11
3305 3191
3306 3192
3307 3193 listMeteors1.append(meteorAux)
3308 3194
3309 3195 return listMeteors1
3310 3196
3311 3197 #Exponential Function
3312 3198
3313 3199 def __exponential_function(self, x, a, tau):
3314 3200 y = a*numpy.exp(-x/tau)
3315 3201 return y
3316 3202
3317 3203 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
3318 3204
3319 3205 pairslist1 = list(pairslist)
3320 3206 pairslist1.append((0,1))
3321 3207 pairslist1.append((3,4))
3322 3208 numPairs = len(pairslist1)
3323 3209 #Time Lag
3324 3210 timeLag = 45*10**-3
3325 3211 c = 3e8
3326 3212 lag = numpy.ceil(timeLag/timeInterval)
3327 3213 freq = 30e6
3328 3214
3329 3215 listMeteors1 = []
3330 3216
3331 3217 for i in range(len(listMeteors)):
3332 3218 meteorAux = listMeteors[i]
3333 3219 if meteorAux[-1] == 0:
3334 3220 mStart = listMeteors[i][1]
3335 3221 mPeak = listMeteors[i][2]
3336 3222 mLag = mPeak - mStart + lag
3337 3223
3338 3224 #get the volt data between the start and end times of the meteor
3339 3225 meteorVolts = listVolts[i]
3340 3226 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3341 3227
3342 3228 #Get CCF
3343 3229 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
3344 3230
3345 3231 #Method 2
3346 3232 slopes = numpy.zeros(numPairs)
3347 3233 time = numpy.array([-2,-1,1,2])*timeInterval
3348 3234 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
3349 3235
3350 3236 #Correct phases
3351 3237 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
3352 3238 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
3353 3239
3354 3240 if indDer[0].shape[0] > 0:
3355 3241 for i in range(indDer[0].shape[0]):
3356 3242 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
3357 3243 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
3358 3244
3359 3245 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
3360 3246 for j in range(numPairs):
3361 3247 fit = stats.linregress(time, angAllCCF[j,:])
3362 3248 slopes[j] = fit[0]
3363 3249
3364 3250 #Remove Outlier
3365 3251 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3366 3252 # slopes = numpy.delete(slopes,indOut)
3367 3253 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3368 3254 # slopes = numpy.delete(slopes,indOut)
3369 3255
3370 3256 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
3371 3257 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
3372 3258 meteorAux[-2] = radialError
3373 3259 meteorAux[-3] = radialVelocity
3374 3260
3375 3261 #Setting Error
3376 3262 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
3377 3263 if numpy.abs(radialVelocity) > 200:
3378 3264 meteorAux[-1] = 15
3379 3265 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
3380 3266 elif radialError > radialStdThresh:
3381 3267 meteorAux[-1] = 12
3382 3268
3383 3269 listMeteors1.append(meteorAux)
3384 3270 return listMeteors1
3385 3271
3386 3272 def __setNewArrays(self, listMeteors, date, heiRang):
3387 3273
3388 3274 #New arrays
3389 3275 arrayMeteors = numpy.array(listMeteors)
3390 3276 arrayParameters = numpy.zeros((len(listMeteors), 13))
3391 3277
3392 3278 #Date inclusion
3393 3279 # date = re.findall(r'\((.*?)\)', date)
3394 3280 # date = date[0].split(',')
3395 3281 # date = map(int, date)
3396 3282 #
3397 3283 # if len(date)<6:
3398 3284 # date.append(0)
3399 3285 #
3400 3286 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
3401 3287 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
3402 3288 arrayDate = numpy.tile(date, (len(listMeteors)))
3403 3289
3404 3290 #Meteor array
3405 3291 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
3406 3292 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
3407 3293
3408 3294 #Parameters Array
3409 3295 arrayParameters[:,0] = arrayDate #Date
3410 3296 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
3411 3297 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
3412 3298 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
3413 3299 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
3414 3300
3415 3301
3416 3302 return arrayParameters
3417 3303
3418 3304 class CorrectSMPhases(Operation):
3419 3305
3420 3306 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
3421 3307
3422 3308 arrayParameters = dataOut.data_param
3423 3309 pairsList = []
3424 3310 pairx = (0,1)
3425 3311 pairy = (2,3)
3426 3312 pairsList.append(pairx)
3427 3313 pairsList.append(pairy)
3428 3314 jph = numpy.zeros(4)
3429 3315
3430 3316 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
3431 3317 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
3432 3318 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
3433 3319
3434 3320 meteorOps = SMOperations()
3435 3321 if channelPositions is None:
3436 3322 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3437 3323 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3438 3324
3439 3325 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3440 3326 h = (hmin,hmax)
3441 3327
3442 3328 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
3443 3329
3444 3330 dataOut.data_param = arrayParameters
3445 3331 return
3446 3332
3447 3333 class SMPhaseCalibration(Operation):
3448 3334
3449 3335 __buffer = None
3450 3336
3451 3337 __initime = None
3452 3338
3453 3339 __dataReady = False
3454 3340
3455 3341 __isConfig = False
3456 3342
3457 3343 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
3458 3344
3459 3345 dataTime = currentTime + paramInterval
3460 3346 deltaTime = dataTime - initTime
3461 3347
3462 3348 if deltaTime >= outputInterval or deltaTime < 0:
3463 3349 return True
3464 3350
3465 3351 return False
3466 3352
3467 3353 def __getGammas(self, pairs, d, phases):
3468 3354 gammas = numpy.zeros(2)
3469 3355
3470 3356 for i in range(len(pairs)):
3471 3357
3472 3358 pairi = pairs[i]
3473 3359
3474 3360 phip3 = phases[:,pairi[0]]
3475 3361 d3 = d[pairi[0]]
3476 3362 phip2 = phases[:,pairi[1]]
3477 3363 d2 = d[pairi[1]]
3478 3364 #Calculating gamma
3479 3365 # jdcos = alp1/(k*d1)
3480 3366 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
3481 3367 jgamma = -phip2*d3/d2 - phip3
3482 3368 jgamma = numpy.angle(numpy.exp(1j*jgamma))
3483 3369 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
3484 3370 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
3485 3371
3486 3372 #Revised distribution
3487 3373 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
3488 3374
3489 3375 #Histogram
3490 3376 nBins = 64
3491 3377 rmin = -0.5*numpy.pi
3492 3378 rmax = 0.5*numpy.pi
3493 3379 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
3494 3380
3495 3381 meteorsY = phaseHisto[0]
3496 3382 phasesX = phaseHisto[1][:-1]
3497 3383 width = phasesX[1] - phasesX[0]
3498 3384 phasesX += width/2
3499 3385
3500 3386 #Gaussian aproximation
3501 3387 bpeak = meteorsY.argmax()
3502 3388 peak = meteorsY.max()
3503 3389 jmin = bpeak - 5
3504 3390 jmax = bpeak + 5 + 1
3505 3391
3506 3392 if jmin<0:
3507 3393 jmin = 0
3508 3394 jmax = 6
3509 3395 elif jmax > meteorsY.size:
3510 3396 jmin = meteorsY.size - 6
3511 3397 jmax = meteorsY.size
3512 3398
3513 3399 x0 = numpy.array([peak,bpeak,50])
3514 3400 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
3515 3401
3516 3402 #Gammas
3517 3403 gammas[i] = coeff[0][1]
3518 3404
3519 3405 return gammas
3520 3406
3521 3407 def __residualFunction(self, coeffs, y, t):
3522 3408
3523 3409 return y - self.__gauss_function(t, coeffs)
3524 3410
3525 3411 def __gauss_function(self, t, coeffs):
3526 3412
3527 3413 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
3528 3414
3529 3415 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
3530 3416 meteorOps = SMOperations()
3531 3417 nchan = 4
3532 3418 pairx = pairsList[0] #x es 0
3533 3419 pairy = pairsList[1] #y es 1
3534 3420 center_xangle = 0
3535 3421 center_yangle = 0
3536 3422 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
3537 3423 ntimes = len(range_angle)
3538 3424
3539 3425 nstepsx = 20
3540 3426 nstepsy = 20
3541 3427
3542 3428 for iz in range(ntimes):
3543 3429 min_xangle = -range_angle[iz]/2 + center_xangle
3544 3430 max_xangle = range_angle[iz]/2 + center_xangle
3545 3431 min_yangle = -range_angle[iz]/2 + center_yangle
3546 3432 max_yangle = range_angle[iz]/2 + center_yangle
3547 3433
3548 3434 inc_x = (max_xangle-min_xangle)/nstepsx
3549 3435 inc_y = (max_yangle-min_yangle)/nstepsy
3550 3436
3551 3437 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
3552 3438 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
3553 3439 penalty = numpy.zeros((nstepsx,nstepsy))
3554 3440 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
3555 3441 jph = numpy.zeros(nchan)
3556 3442
3557 3443 # Iterations looking for the offset
3558 3444 for iy in range(int(nstepsy)):
3559 3445 for ix in range(int(nstepsx)):
3560 3446 d3 = d[pairsList[1][0]]
3561 3447 d2 = d[pairsList[1][1]]
3562 3448 d5 = d[pairsList[0][0]]
3563 3449 d4 = d[pairsList[0][1]]
3564 3450
3565 3451 alp2 = alpha_y[iy] #gamma 1
3566 3452 alp4 = alpha_x[ix] #gamma 0
3567 3453
3568 3454 alp3 = -alp2*d3/d2 - gammas[1]
3569 3455 alp5 = -alp4*d5/d4 - gammas[0]
3570 3456 # jph[pairy[1]] = alpha_y[iy]
3571 3457 # jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
3572 3458
3573 3459 # jph[pairx[1]] = alpha_x[ix]
3574 3460 # jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
3575 3461 jph[pairsList[0][1]] = alp4
3576 3462 jph[pairsList[0][0]] = alp5
3577 3463 jph[pairsList[1][0]] = alp3
3578 3464 jph[pairsList[1][1]] = alp2
3579 3465 jph_array[:,ix,iy] = jph
3580 3466 # d = [2.0,2.5,2.5,2.0]
3581 3467 #falta chequear si va a leer bien los meteoros
3582 3468 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
3583 3469 error = meteorsArray1[:,-1]
3584 3470 ind1 = numpy.where(error==0)[0]
3585 3471 penalty[ix,iy] = ind1.size
3586 3472
3587 3473 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
3588 3474 phOffset = jph_array[:,i,j]
3589 3475
3590 3476 center_xangle = phOffset[pairx[1]]
3591 3477 center_yangle = phOffset[pairy[1]]
3592 3478
3593 3479 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
3594 3480 phOffset = phOffset*180/numpy.pi
3595 3481 return phOffset
3596 3482
3597 3483
3598 3484 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
3599 3485
3600 3486 dataOut.flagNoData = True
3601 3487 self.__dataReady = False
3602 3488 dataOut.outputInterval = nHours*3600
3603 3489
3604 3490 if self.__isConfig == False:
3605 3491 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
3606 3492 #Get Initial LTC time
3607 3493 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3608 3494 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
3609 3495
3610 3496 self.__isConfig = True
3611 3497
3612 3498 if self.__buffer is None:
3613 3499 self.__buffer = dataOut.data_param.copy()
3614 3500
3615 3501 else:
3616 3502 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
3617 3503
3618 3504 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
3619 3505
3620 3506 if self.__dataReady:
3621 3507 dataOut.utctimeInit = self.__initime
3622 3508 self.__initime += dataOut.outputInterval #to erase time offset
3623 3509
3624 3510 freq = dataOut.frequency
3625 3511 c = dataOut.C #m/s
3626 3512 lamb = c/freq
3627 3513 k = 2*numpy.pi/lamb
3628 3514 azimuth = 0
3629 3515 h = (hmin, hmax)
3630 3516 # pairs = ((0,1),(2,3)) #Estrella
3631 3517 # pairs = ((1,0),(2,3)) #T
3632 3518
3633 3519 if channelPositions is None:
3634 3520 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3635 3521 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3636 3522 meteorOps = SMOperations()
3637 3523 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3638 3524
3639 3525 #Checking correct order of pairs
3640 3526 pairs = []
3641 3527 if distances[1] > distances[0]:
3642 3528 pairs.append((1,0))
3643 3529 else:
3644 3530 pairs.append((0,1))
3645 3531
3646 3532 if distances[3] > distances[2]:
3647 3533 pairs.append((3,2))
3648 3534 else:
3649 3535 pairs.append((2,3))
3650 3536 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
3651 3537
3652 3538 meteorsArray = self.__buffer
3653 3539 error = meteorsArray[:,-1]
3654 3540 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
3655 3541 ind1 = numpy.where(boolError)[0]
3656 3542 meteorsArray = meteorsArray[ind1,:]
3657 3543 meteorsArray[:,-1] = 0
3658 3544 phases = meteorsArray[:,8:12]
3659 3545
3660 3546 #Calculate Gammas
3661 3547 gammas = self.__getGammas(pairs, distances, phases)
3662 3548 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
3663 3549 #Calculate Phases
3664 3550 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
3665 3551 phasesOff = phasesOff.reshape((1,phasesOff.size))
3666 3552 dataOut.data_output = -phasesOff
3667 3553 dataOut.flagNoData = False
3668 3554 self.__buffer = None
3669 3555
3670 3556
3671 3557 return
3672 3558
3673 3559 class SMOperations():
3674 3560
3675 3561 def __init__(self):
3676 3562
3677 3563 return
3678 3564
3679 3565 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
3680 3566
3681 3567 arrayParameters = arrayParameters0.copy()
3682 3568 hmin = h[0]
3683 3569 hmax = h[1]
3684 3570
3685 3571 #Calculate AOA (Error N 3, 4)
3686 3572 #JONES ET AL. 1998
3687 3573 AOAthresh = numpy.pi/8
3688 3574 error = arrayParameters[:,-1]
3689 3575 phases = -arrayParameters[:,8:12] + jph
3690 3576 # phases = numpy.unwrap(phases)
3691 3577 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
3692 3578
3693 3579 #Calculate Heights (Error N 13 and 14)
3694 3580 error = arrayParameters[:,-1]
3695 3581 Ranges = arrayParameters[:,1]
3696 3582 zenith = arrayParameters[:,4]
3697 3583 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
3698 3584
3699 3585 #----------------------- Get Final data ------------------------------------
3700 3586 # error = arrayParameters[:,-1]
3701 3587 # ind1 = numpy.where(error==0)[0]
3702 3588 # arrayParameters = arrayParameters[ind1,:]
3703 3589
3704 3590 return arrayParameters
3705 3591
3706 3592 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
3707 3593
3708 3594 arrayAOA = numpy.zeros((phases.shape[0],3))
3709 3595 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
3710 3596
3711 3597 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3712 3598 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3713 3599 arrayAOA[:,2] = cosDirError
3714 3600
3715 3601 azimuthAngle = arrayAOA[:,0]
3716 3602 zenithAngle = arrayAOA[:,1]
3717 3603
3718 3604 #Setting Error
3719 3605 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
3720 3606 error[indError] = 0
3721 3607 #Number 3: AOA not fesible
3722 3608 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3723 3609 error[indInvalid] = 3
3724 3610 #Number 4: Large difference in AOAs obtained from different antenna baselines
3725 3611 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3726 3612 error[indInvalid] = 4
3727 3613 return arrayAOA, error
3728 3614
3729 3615 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
3730 3616
3731 3617 #Initializing some variables
3732 3618 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3733 3619 ang_aux = ang_aux.reshape(1,ang_aux.size)
3734 3620
3735 3621 cosdir = numpy.zeros((arrayPhase.shape[0],2))
3736 3622 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3737 3623
3738 3624
3739 3625 for i in range(2):
3740 3626 ph0 = arrayPhase[:,pairsList[i][0]]
3741 3627 ph1 = arrayPhase[:,pairsList[i][1]]
3742 3628 d0 = distances[pairsList[i][0]]
3743 3629 d1 = distances[pairsList[i][1]]
3744 3630
3745 3631 ph0_aux = ph0 + ph1
3746 3632 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
3747 3633 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
3748 3634 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
3749 3635 #First Estimation
3750 3636 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
3751 3637
3752 3638 #Most-Accurate Second Estimation
3753 3639 phi1_aux = ph0 - ph1
3754 3640 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3755 3641 #Direction Cosine 1
3756 3642 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
3757 3643
3758 3644 #Searching the correct Direction Cosine
3759 3645 cosdir0_aux = cosdir0[:,i]
3760 3646 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3761 3647 #Minimum Distance
3762 3648 cosDiff = (cosdir1 - cosdir0_aux)**2
3763 3649 indcos = cosDiff.argmin(axis = 1)
3764 3650 #Saving Value obtained
3765 3651 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3766 3652
3767 3653 return cosdir0, cosdir
3768 3654
3769 3655 def __calculateAOA(self, cosdir, azimuth):
3770 3656 cosdirX = cosdir[:,0]
3771 3657 cosdirY = cosdir[:,1]
3772 3658
3773 3659 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3774 3660 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
3775 3661 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3776 3662
3777 3663 return angles
3778 3664
3779 3665 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3780 3666
3781 3667 Ramb = 375 #Ramb = c/(2*PRF)
3782 3668 Re = 6371 #Earth Radius
3783 3669 heights = numpy.zeros(Ranges.shape)
3784 3670
3785 3671 R_aux = numpy.array([0,1,2])*Ramb
3786 3672 R_aux = R_aux.reshape(1,R_aux.size)
3787 3673
3788 3674 Ranges = Ranges.reshape(Ranges.size,1)
3789 3675
3790 3676 Ri = Ranges + R_aux
3791 3677 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3792 3678
3793 3679 #Check if there is a height between 70 and 110 km
3794 3680 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3795 3681 ind_h = numpy.where(h_bool == 1)[0]
3796 3682
3797 3683 hCorr = hi[ind_h, :]
3798 3684 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3799 3685
3800 3686 hCorr = hi[ind_hCorr][:len(ind_h)]
3801 3687 heights[ind_h] = hCorr
3802 3688
3803 3689 #Setting Error
3804 3690 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3805 3691 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3806 3692 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
3807 3693 error[indError] = 0
3808 3694 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3809 3695 error[indInvalid2] = 14
3810 3696 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3811 3697 error[indInvalid1] = 13
3812 3698
3813 3699 return heights, error
3814 3700
3815 3701 def getPhasePairs(self, channelPositions):
3816 3702 chanPos = numpy.array(channelPositions)
3817 3703 listOper = list(itertools.combinations(list(range(5)),2))
3818 3704
3819 3705 distances = numpy.zeros(4)
3820 3706 axisX = []
3821 3707 axisY = []
3822 3708 distX = numpy.zeros(3)
3823 3709 distY = numpy.zeros(3)
3824 3710 ix = 0
3825 3711 iy = 0
3826 3712
3827 3713 pairX = numpy.zeros((2,2))
3828 3714 pairY = numpy.zeros((2,2))
3829 3715
3830 3716 for i in range(len(listOper)):
3831 3717 pairi = listOper[i]
3832 3718
3833 3719 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
3834 3720
3835 3721 if posDif[0] == 0:
3836 3722 axisY.append(pairi)
3837 3723 distY[iy] = posDif[1]
3838 3724 iy += 1
3839 3725 elif posDif[1] == 0:
3840 3726 axisX.append(pairi)
3841 3727 distX[ix] = posDif[0]
3842 3728 ix += 1
3843 3729
3844 3730 for i in range(2):
3845 3731 if i==0:
3846 3732 dist0 = distX
3847 3733 axis0 = axisX
3848 3734 else:
3849 3735 dist0 = distY
3850 3736 axis0 = axisY
3851 3737
3852 3738 side = numpy.argsort(dist0)[:-1]
3853 3739 axis0 = numpy.array(axis0)[side,:]
3854 3740 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
3855 3741 axis1 = numpy.unique(numpy.reshape(axis0,4))
3856 3742 side = axis1[axis1 != chanC]
3857 3743 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
3858 3744 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
3859 3745 if diff1<0:
3860 3746 chan2 = side[0]
3861 3747 d2 = numpy.abs(diff1)
3862 3748 chan1 = side[1]
3863 3749 d1 = numpy.abs(diff2)
3864 3750 else:
3865 3751 chan2 = side[1]
3866 3752 d2 = numpy.abs(diff2)
3867 3753 chan1 = side[0]
3868 3754 d1 = numpy.abs(diff1)
3869 3755
3870 3756 if i==0:
3871 3757 chanCX = chanC
3872 3758 chan1X = chan1
3873 3759 chan2X = chan2
3874 3760 distances[0:2] = numpy.array([d1,d2])
3875 3761 else:
3876 3762 chanCY = chanC
3877 3763 chan1Y = chan1
3878 3764 chan2Y = chan2
3879 3765 distances[2:4] = numpy.array([d1,d2])
3880 3766 # axisXsides = numpy.reshape(axisX[ix,:],4)
3881 3767 #
3882 3768 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
3883 3769 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
3884 3770 #
3885 3771 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
3886 3772 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
3887 3773 # channel25X = int(pairX[0,ind25X])
3888 3774 # channel20X = int(pairX[1,ind20X])
3889 3775 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
3890 3776 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
3891 3777 # channel25Y = int(pairY[0,ind25Y])
3892 3778 # channel20Y = int(pairY[1,ind20Y])
3893 3779
3894 3780 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
3895 3781 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
3896 3782
3897 3783 return pairslist, distances
3898 3784 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
3899 3785 #
3900 3786 # arrayAOA = numpy.zeros((phases.shape[0],3))
3901 3787 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
3902 3788 #
3903 3789 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3904 3790 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3905 3791 # arrayAOA[:,2] = cosDirError
3906 3792 #
3907 3793 # azimuthAngle = arrayAOA[:,0]
3908 3794 # zenithAngle = arrayAOA[:,1]
3909 3795 #
3910 3796 # #Setting Error
3911 3797 # #Number 3: AOA not fesible
3912 3798 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3913 3799 # error[indInvalid] = 3
3914 3800 # #Number 4: Large difference in AOAs obtained from different antenna baselines
3915 3801 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3916 3802 # error[indInvalid] = 4
3917 3803 # return arrayAOA, error
3918 3804 #
3919 3805 # def __getDirectionCosines(self, arrayPhase, pairsList):
3920 3806 #
3921 3807 # #Initializing some variables
3922 3808 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3923 3809 # ang_aux = ang_aux.reshape(1,ang_aux.size)
3924 3810 #
3925 3811 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
3926 3812 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3927 3813 #
3928 3814 #
3929 3815 # for i in range(2):
3930 3816 # #First Estimation
3931 3817 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
3932 3818 # #Dealias
3933 3819 # indcsi = numpy.where(phi0_aux > numpy.pi)
3934 3820 # phi0_aux[indcsi] -= 2*numpy.pi
3935 3821 # indcsi = numpy.where(phi0_aux < -numpy.pi)
3936 3822 # phi0_aux[indcsi] += 2*numpy.pi
3937 3823 # #Direction Cosine 0
3938 3824 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
3939 3825 #
3940 3826 # #Most-Accurate Second Estimation
3941 3827 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
3942 3828 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3943 3829 # #Direction Cosine 1
3944 3830 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
3945 3831 #
3946 3832 # #Searching the correct Direction Cosine
3947 3833 # cosdir0_aux = cosdir0[:,i]
3948 3834 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3949 3835 # #Minimum Distance
3950 3836 # cosDiff = (cosdir1 - cosdir0_aux)**2
3951 3837 # indcos = cosDiff.argmin(axis = 1)
3952 3838 # #Saving Value obtained
3953 3839 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3954 3840 #
3955 3841 # return cosdir0, cosdir
3956 3842 #
3957 3843 # def __calculateAOA(self, cosdir, azimuth):
3958 3844 # cosdirX = cosdir[:,0]
3959 3845 # cosdirY = cosdir[:,1]
3960 3846 #
3961 3847 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3962 3848 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
3963 3849 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3964 3850 #
3965 3851 # return angles
3966 3852 #
3967 3853 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3968 3854 #
3969 3855 # Ramb = 375 #Ramb = c/(2*PRF)
3970 3856 # Re = 6371 #Earth Radius
3971 3857 # heights = numpy.zeros(Ranges.shape)
3972 3858 #
3973 3859 # R_aux = numpy.array([0,1,2])*Ramb
3974 3860 # R_aux = R_aux.reshape(1,R_aux.size)
3975 3861 #
3976 3862 # Ranges = Ranges.reshape(Ranges.size,1)
3977 3863 #
3978 3864 # Ri = Ranges + R_aux
3979 3865 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3980 3866 #
3981 3867 # #Check if there is a height between 70 and 110 km
3982 3868 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3983 3869 # ind_h = numpy.where(h_bool == 1)[0]
3984 3870 #
3985 3871 # hCorr = hi[ind_h, :]
3986 3872 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3987 3873 #
3988 3874 # hCorr = hi[ind_hCorr]
3989 3875 # heights[ind_h] = hCorr
3990 3876 #
3991 3877 # #Setting Error
3992 3878 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3993 3879 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3994 3880 #
3995 3881 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3996 3882 # error[indInvalid2] = 14
3997 3883 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3998 3884 # error[indInvalid1] = 13
3999 3885 #
4000 3886 # return heights, error
General Comments 0
You need to be logged in to leave comments. Login now