|
This diff has been collapsed as it changes many lines, (2356 lines changed)
Show them
Hide them
|
|
|
@@
-1,6
+1,7
|
|
|
1
|
import numpy
|
|
1
|
import numpy
|
|
|
2
|
import math
|
|
2
|
import math
|
|
|
3
|
from scipy import optimize, interpolate, signal, stats, ndimage
|
|
3
|
from scipy import optimize, interpolate, signal, stats, ndimage
|
|
|
|
|
|
4
|
from scipy.stats import norm
|
|
|
4
|
import scipy
|
|
5
|
import scipy
|
|
|
5
|
import re
|
|
6
|
import re
|
|
|
6
|
import datetime
|
|
7
|
import datetime
|
|
|
@@
-13,20
+14,20
from multiprocessing.pool import ThreadPool
|
|
|
13
|
import time
|
|
14
|
import time
|
|
|
14
|
|
|
15
|
|
|
|
15
|
from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
|
|
16
|
from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
|
|
|
16
|
from .jroproc_base import ProcessingUnit, Operation, MPDecorator
|
|
17
|
from .jroproc_base import ProcessingUnit, Operation #, MPDecorator
|
|
|
17
|
from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
|
|
18
|
from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
|
|
|
18
|
from scipy import asarray as ar,exp
|
|
19
|
from scipy import asarray as ar,exp
|
|
|
19
|
from scipy.optimize import curve_fit
|
|
20
|
from scipy.optimize import curve_fit
|
|
|
20
|
from schainpy.utils import log
|
|
21
|
#from schainpy.utils import log
|
|
|
21
|
import warnings
|
|
22
|
import warnings
|
|
|
22
|
from numpy import NaN
|
|
23
|
from numpy import NaN
|
|
|
23
|
from scipy.optimize.optimize import OptimizeWarning
|
|
24
|
from scipy.optimize.optimize import OptimizeWarning
|
|
|
24
|
warnings.filterwarnings('ignore')
|
|
25
|
warnings.filterwarnings('ignore')
|
|
|
25
|
|
|
26
|
|
|
|
26
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
27
|
|
|
27
|
|
|
|
28
|
SPEED_OF_LIGHT = 299792458
|
|
28
|
SPEED_OF_LIGHT = 299792458
|
|
|
29
|
|
|
29
|
|
|
|
|
|
|
30
|
|
|
|
30
|
'''solving pickling issue'''
|
|
31
|
'''solving pickling issue'''
|
|
|
31
|
|
|
32
|
|
|
|
32
|
def _pickle_method(method):
|
|
33
|
def _pickle_method(method):
|
|
|
@@
-45,7
+46,7
def _unpickle_method(func_name, obj, cls):
|
|
|
45
|
break
|
|
46
|
break
|
|
|
46
|
return func.__get__(obj, cls)
|
|
47
|
return func.__get__(obj, cls)
|
|
|
47
|
|
|
48
|
|
|
|
48
|
|
|
49
|
#@MPDecorator
|
|
|
49
|
class ParametersProc(ProcessingUnit):
|
|
50
|
class ParametersProc(ProcessingUnit):
|
|
|
50
|
|
|
51
|
|
|
|
51
|
METHODS = {}
|
|
52
|
METHODS = {}
|
|
|
@@
-77,9
+78,9
class ParametersProc(ProcessingUnit):
|
|
|
77
|
self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
|
|
78
|
self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
|
|
|
78
|
#self.dataOut.nHeights = self.dataIn.nHeights
|
|
79
|
#self.dataOut.nHeights = self.dataIn.nHeights
|
|
|
79
|
#self.dataOut.nChannels = self.dataIn.nChannels
|
|
80
|
#self.dataOut.nChannels = self.dataIn.nChannels
|
|
|
80
|
# self.dataOut.nBaud = self.dataIn.nBaud
|
|
81
|
self.dataOut.nBaud = self.dataIn.nBaud
|
|
|
81
|
# self.dataOut.nCode = self.dataIn.nCode
|
|
82
|
self.dataOut.nCode = self.dataIn.nCode
|
|
|
82
|
# self.dataOut.code = self.dataIn.code
|
|
83
|
self.dataOut.code = self.dataIn.code
|
|
|
83
|
#self.dataOut.nProfiles = self.dataOut.nFFTPoints
|
|
84
|
#self.dataOut.nProfiles = self.dataOut.nFFTPoints
|
|
|
84
|
self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
|
|
85
|
self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
|
|
|
85
|
# self.dataOut.utctime = self.firstdatatime
|
|
86
|
# self.dataOut.utctime = self.firstdatatime
|
|
|
@@
-88,10
+89,10
class ParametersProc(ProcessingUnit):
|
|
|
88
|
self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
|
|
89
|
self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
|
|
|
89
|
self.dataOut.nCohInt = self.dataIn.nCohInt
|
|
90
|
self.dataOut.nCohInt = self.dataIn.nCohInt
|
|
|
90
|
#self.dataOut.nIncohInt = 1
|
|
91
|
#self.dataOut.nIncohInt = 1
|
|
|
91
|
# self.dataOut.ippSeconds = self.dataIn.ippSeconds
|
|
92
|
self.dataOut.ippSeconds = self.dataIn.ippSeconds
|
|
|
92
|
# self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
|
|
93
|
# self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
|
|
|
93
|
self.dataOut.timeInterval1 = self.dataIn.timeInterval
|
|
94
|
self.dataOut.timeInterval1 = self.dataIn.timeInterval
|
|
|
94
|
self.dataOut.heightList = self.dataIn.heightList
|
|
95
|
self.dataOut.heightList = self.dataIn.heightList #getHeiRange()
|
|
|
95
|
self.dataOut.frequency = self.dataIn.frequency
|
|
96
|
self.dataOut.frequency = self.dataIn.frequency
|
|
|
96
|
#self.dataOut.noise = self.dataIn.noise
|
|
97
|
#self.dataOut.noise = self.dataIn.noise
|
|
|
97
|
|
|
98
|
|
|
|
@@
-108,27
+109,13
class ParametersProc(ProcessingUnit):
|
|
|
108
|
self.dataOut.flagNoData = False
|
|
109
|
self.dataOut.flagNoData = False
|
|
|
109
|
self.dataOut.utctimeInit = self.dataIn.utctime
|
|
110
|
self.dataOut.utctimeInit = self.dataIn.utctime
|
|
|
110
|
self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
|
|
111
|
self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
|
|
|
111
|
if hasattr(self.dataIn, 'dataPP_POW'):
|
|
|
|
|
|
112
|
self.dataOut.dataPP_POW = self.dataIn.dataPP_POW
|
|
|
|
|
|
113
|
|
|
|
|
|
|
114
|
if hasattr(self.dataIn, 'dataPP_POWER'):
|
|
|
|
|
|
115
|
self.dataOut.dataPP_POWER = self.dataIn.dataPP_POWER
|
|
|
|
|
|
116
|
|
|
|
|
|
|
117
|
if hasattr(self.dataIn, 'dataPP_DOP'):
|
|
|
|
|
|
118
|
self.dataOut.dataPP_DOP = self.dataIn.dataPP_DOP
|
|
|
|
|
|
119
|
|
|
|
|
|
|
120
|
if hasattr(self.dataIn, 'dataPP_SNR'):
|
|
|
|
|
|
121
|
self.dataOut.dataPP_SNR = self.dataIn.dataPP_SNR
|
|
|
|
|
|
122
|
|
|
|
|
|
|
123
|
if hasattr(self.dataIn, 'dataPP_WIDTH'):
|
|
|
|
|
|
124
|
self.dataOut.dataPP_WIDTH = self.dataIn.dataPP_WIDTH
|
|
|
|
|
|
125
|
return
|
|
112
|
return
|
|
|
126
|
|
|
113
|
|
|
|
127
|
#---------------------- Spectra Data ---------------------------
|
|
114
|
#---------------------- Spectra Data ---------------------------
|
|
|
128
|
|
|
115
|
|
|
|
129
|
if self.dataIn.type == "Spectra":
|
|
116
|
if self.dataIn.type == "Spectra":
|
|
|
130
|
|
|
117
|
|
|
|
131
|
self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
|
|
118
|
self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
|
|
|
132
|
self.dataOut.data_spc = self.dataIn.data_spc
|
|
119
|
self.dataOut.data_spc = self.dataIn.data_spc
|
|
|
133
|
self.dataOut.data_cspc = self.dataIn.data_cspc
|
|
120
|
self.dataOut.data_cspc = self.dataIn.data_cspc
|
|
|
134
|
self.dataOut.nProfiles = self.dataIn.nProfiles
|
|
121
|
self.dataOut.nProfiles = self.dataIn.nProfiles
|
|
|
@@
-142,6
+129,7
class ParametersProc(ProcessingUnit):
|
|
|
142
|
self.dataOut.pairsList = self.dataIn.pairsList
|
|
129
|
self.dataOut.pairsList = self.dataIn.pairsList
|
|
|
143
|
self.dataOut.groupList = self.dataIn.pairsList
|
|
130
|
self.dataOut.groupList = self.dataIn.pairsList
|
|
|
144
|
self.dataOut.flagNoData = False
|
|
131
|
self.dataOut.flagNoData = False
|
|
|
|
|
|
132
|
self.dataOut.spcacum = None
|
|
|
145
|
|
|
133
|
|
|
|
146
|
if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
|
|
134
|
if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
|
|
|
147
|
self.dataOut.ChanDist = self.dataIn.ChanDist
|
|
135
|
self.dataOut.ChanDist = self.dataIn.ChanDist
|
|
|
@@
-173,7
+161,7
class ParametersProc(ProcessingUnit):
|
|
|
173
|
|
|
161
|
|
|
|
174
|
self.dataOut.abscissaList = self.dataIn.lagRange
|
|
162
|
self.dataOut.abscissaList = self.dataIn.lagRange
|
|
|
175
|
self.dataOut.noise = self.dataIn.noise
|
|
163
|
self.dataOut.noise = self.dataIn.noise
|
|
|
176
|
self.dataOut.data_snr = self.dataIn.SNR
|
|
164
|
self.dataOut.data_SNR = self.dataIn.SNR
|
|
|
177
|
self.dataOut.flagNoData = False
|
|
165
|
self.dataOut.flagNoData = False
|
|
|
178
|
self.dataOut.nAvg = self.dataIn.nAvg
|
|
166
|
self.dataOut.nAvg = self.dataIn.nAvg
|
|
|
179
|
|
|
167
|
|
|
|
@@
-198,11
+186,13
def target(tups):
|
|
|
198
|
|
|
186
|
|
|
|
199
|
return obj.FitGau(args)
|
|
187
|
return obj.FitGau(args)
|
|
|
200
|
|
|
188
|
|
|
|
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
|
|
|
|
|
|
204
|
|
|
189
|
|
|
|
205
|
ClutterWidth : Width to look for the clutter peak
|
|
190
|
class SpectralFilters(Operation):
|
|
|
|
|
|
191
|
|
|
|
|
|
|
192
|
'''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
|
|
|
|
|
|
193
|
|
|
|
|
|
|
194
|
LimitR : It is the limit in m/s of Rainfall
|
|
|
|
|
|
195
|
LimitW : It is the limit in m/s for Winds
|
|
|
206
|
|
|
196
|
|
|
|
207
|
Input:
|
|
197
|
Input:
|
|
|
208
|
|
|
198
|
|
|
|
@@
-212,111
+202,91
class RemoveWideGC(Operation):
|
|
|
212
|
Affected:
|
|
202
|
Affected:
|
|
|
213
|
|
|
203
|
|
|
|
214
|
self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
|
|
204
|
self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
|
|
|
|
|
|
205
|
self.dataOut.spcparam_range : Used in SpcParamPlot
|
|
|
|
|
|
206
|
self.dataOut.SPCparam : Used in PrecipitationProc
|
|
|
|
|
|
207
|
|
|
|
215
|
|
|
208
|
|
|
|
216
|
Written by D. Scipión 25.02.2021
|
|
|
|
|
|
217
|
'''
|
|
209
|
'''
|
|
|
|
|
|
210
|
|
|
|
218
|
def __init__(self):
|
|
211
|
def __init__(self):
|
|
|
219
|
Operation.__init__(self)
|
|
212
|
Operation.__init__(self)
|
|
|
220
|
self.i=0
|
|
213
|
self.i=0
|
|
|
221
|
self.ich = 0
|
|
|
|
|
|
222
|
self.ir = 0
|
|
|
|
|
|
223
|
|
|
214
|
|
|
|
224
|
def run(self, dataOut, ClutterWidth=2.5):
|
|
215
|
def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
|
|
|
225
|
# print ('Entering RemoveWideGC ... ')
|
|
216
|
|
|
|
|
|
|
217
|
|
|
|
|
|
|
218
|
#Limite de vientos
|
|
|
|
|
|
219
|
LimitR = PositiveLimit
|
|
|
|
|
|
220
|
LimitN = NegativeLimit
|
|
|
226
|
|
|
221
|
|
|
|
227
|
self.spc = dataOut.data_pre[0].copy()
|
|
222
|
self.spc = dataOut.data_pre[0].copy()
|
|
|
228
|
self.spc_out = dataOut.data_pre[0].copy()
|
|
223
|
self.cspc = dataOut.data_pre[1].copy()
|
|
|
229
|
self.Num_Chn = self.spc.shape[0]
|
|
224
|
|
|
|
230
|
self.Num_Hei = self.spc.shape[2]
|
|
225
|
self.Num_Hei = self.spc.shape[2]
|
|
|
231
|
VelRange = dataOut.spc_range[2][:-1]
|
|
226
|
self.Num_Bin = self.spc.shape[1]
|
|
|
232
|
dv = VelRange[1]-VelRange[0]
|
|
227
|
self.Num_Chn = self.spc.shape[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
|
|
|
228
|
|
|
|
255
|
valleyindex = j2index[numpy.where(junk4>1)]
|
|
229
|
VelRange = dataOut.spc_range[2]
|
|
|
256
|
peakindex = j1index[numpy.where(junk3>1)]
|
|
230
|
TimeRange = dataOut.spc_range[1]
|
|
|
|
|
|
231
|
FrecRange = dataOut.spc_range[0]
|
|
|
257
|
|
|
232
|
|
|
|
258
|
isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
|
|
233
|
Vmax= 2*numpy.max(dataOut.spc_range[2])
|
|
|
259
|
if numpy.size(isvalid) == 0 :
|
|
234
|
Tmax= 2*numpy.max(dataOut.spc_range[1])
|
|
|
260
|
continue
|
|
235
|
Fmax= 2*numpy.max(dataOut.spc_range[0])
|
|
|
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
|
|
|
236
|
|
|
|
276
|
# Removing the clutter
|
|
237
|
Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
|
|
|
277
|
interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
|
|
238
|
Breaker1R=numpy.where(VelRange == Breaker1R)
|
|
|
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
|
|
|
239
|
|
|
|
281
|
dataOut.data_pre[0] = self.spc_out
|
|
240
|
Delta = self.Num_Bin/2 - Breaker1R[0]
|
|
|
282
|
#print ('Leaving RemoveWideGC ... ')
|
|
|
|
|
|
283
|
return dataOut
|
|
|
|
|
|
284
|
|
|
241
|
|
|
|
285
|
class SpectralFilters(Operation):
|
|
|
|
|
|
286
|
''' This class allows to replace the novalid values with noise for each channel
|
|
|
|
|
|
287
|
This applies to CLAIRE RADAR
|
|
|
|
|
|
288
|
|
|
242
|
|
|
|
289
|
PositiveLimit : RightLimit of novalid data
|
|
243
|
'''Reacomodando SPCrange'''
|
|
|
290
|
NegativeLimit : LeftLimit of novalid data
|
|
|
|
|
|
291
|
|
|
244
|
|
|
|
292
|
Input:
|
|
245
|
VelRange=numpy.roll(VelRange,-(int(self.Num_Bin/2)) ,axis=0)
|
|
|
293
|
|
|
246
|
|
|
|
294
|
self.dataOut.data_pre : SPC and CSPC
|
|
247
|
VelRange[-(int(self.Num_Bin/2)):]+= Vmax
|
|
|
295
|
self.dataOut.spc_range : To select wind and rainfall velocities
|
|
|
|
|
|
296
|
|
|
248
|
|
|
|
297
|
Affected:
|
|
249
|
FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0)
|
|
|
298
|
|
|
250
|
|
|
|
299
|
self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
|
|
251
|
FrecRange[-(int(self.Num_Bin/2)):]+= Fmax
|
|
|
300
|
|
|
252
|
|
|
|
301
|
Written by D. Scipión 29.01.2021
|
|
253
|
TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0)
|
|
|
302
|
'''
|
|
|
|
|
|
303
|
def __init__(self):
|
|
|
|
|
|
304
|
Operation.__init__(self)
|
|
|
|
|
|
305
|
self.i = 0
|
|
|
|
|
|
306
|
|
|
254
|
|
|
|
307
|
def run(self, dataOut, ):
|
|
255
|
TimeRange[-(int(self.Num_Bin/2)):]+= Tmax
|
|
|
308
|
|
|
256
|
|
|
|
309
|
self.spc = dataOut.data_pre[0].copy()
|
|
257
|
''' ------------------ '''
|
|
|
310
|
self.Num_Chn = self.spc.shape[0]
|
|
258
|
|
|
|
311
|
VelRange = dataOut.spc_range[2]
|
|
259
|
Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
|
|
|
|
|
|
260
|
Breaker2R=numpy.where(VelRange == Breaker2R)
|
|
|
312
|
|
|
261
|
|
|
|
313
|
# novalid corresponds to data within the Negative and PositiveLimit
|
|
|
|
|
|
314
|
|
|
262
|
|
|
|
|
|
|
263
|
SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1)
|
|
|
315
|
|
|
264
|
|
|
|
316
|
# Removing novalid data from the spectra
|
|
265
|
SPCcut = SPCroll.copy()
|
|
|
317
|
for i in range(self.Num_Chn):
|
|
266
|
for i in range(self.Num_Chn):
|
|
|
318
|
self.spc[i,novalid,:] = dataOut.noise[i]
|
|
267
|
|
|
|
319
|
dataOut.data_pre[0] = self.spc
|
|
268
|
SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
|
|
|
|
|
|
269
|
SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
|
|
|
|
|
|
270
|
|
|
|
|
|
|
271
|
SPCcut[i]=SPCcut[i]- dataOut.noise[i]
|
|
|
|
|
|
272
|
SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
|
|
|
|
|
|
273
|
|
|
|
|
|
|
274
|
SPCroll[i]=SPCroll[i]-dataOut.noise[i]
|
|
|
|
|
|
275
|
SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
|
|
|
|
|
|
276
|
|
|
|
|
|
|
277
|
SPC_ch1 = SPCroll
|
|
|
|
|
|
278
|
|
|
|
|
|
|
279
|
SPC_ch2 = SPCcut
|
|
|
|
|
|
280
|
|
|
|
|
|
|
281
|
SPCparam = (SPC_ch1, SPC_ch2, self.spc)
|
|
|
|
|
|
282
|
dataOut.SPCparam = numpy.asarray(SPCparam)
|
|
|
|
|
|
283
|
|
|
|
|
|
|
284
|
|
|
|
|
|
|
285
|
dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
|
|
|
|
|
|
286
|
|
|
|
|
|
|
287
|
dataOut.spcparam_range[2]=VelRange
|
|
|
|
|
|
288
|
dataOut.spcparam_range[1]=TimeRange
|
|
|
|
|
|
289
|
dataOut.spcparam_range[0]=FrecRange
|
|
|
320
|
return dataOut
|
|
290
|
return dataOut
|
|
|
321
|
|
|
291
|
|
|
|
322
|
class GaussianFit(Operation):
|
|
292
|
class GaussianFit(Operation):
|
|
|
@@
-338,163
+308,113
class GaussianFit(Operation):
|
|
|
338
|
self.i=0
|
|
308
|
self.i=0
|
|
|
339
|
|
|
309
|
|
|
|
340
|
|
|
310
|
|
|
|
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
|
|
311
|
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'):
|
|
|
|
|
|
343
|
"""This routine will find a couple of generalized Gaussians to a power spectrum
|
|
312
|
"""This routine will find a couple of generalized Gaussians to a power spectrum
|
|
|
344
|
methods: generalized, squared
|
|
|
|
|
|
345
|
input: spc
|
|
313
|
input: spc
|
|
|
346
|
output:
|
|
314
|
output:
|
|
|
347
|
noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
|
|
315
|
Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
|
|
|
348
|
"""
|
|
316
|
"""
|
|
|
349
|
print ('Entering ',method,' double Gaussian fit')
|
|
317
|
|
|
|
350
|
self.spc = dataOut.data_pre[0].copy()
|
|
318
|
self.spc = dataOut.data_pre[0].copy()
|
|
|
351
|
self.Num_Hei = self.spc.shape[2]
|
|
319
|
self.Num_Hei = self.spc.shape[2]
|
|
|
352
|
self.Num_Bin = self.spc.shape[1]
|
|
320
|
self.Num_Bin = self.spc.shape[1]
|
|
|
353
|
self.Num_Chn = self.spc.shape[0]
|
|
321
|
self.Num_Chn = self.spc.shape[0]
|
|
|
|
|
|
322
|
Vrange = dataOut.abscissaList
|
|
|
|
|
|
323
|
|
|
|
|
|
|
324
|
GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
|
|
|
|
|
|
325
|
SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
|
|
|
|
|
|
326
|
SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
|
|
|
|
|
|
327
|
SPC_ch1[:] = numpy.NaN
|
|
|
|
|
|
328
|
SPC_ch2[:] = numpy.NaN
|
|
|
|
|
|
329
|
|
|
|
354
|
|
|
330
|
|
|
|
355
|
start_time = time.time()
|
|
331
|
start_time = time.time()
|
|
|
356
|
|
|
332
|
|
|
|
|
|
|
333
|
noise_ = dataOut.spc_noise[0].copy()
|
|
|
|
|
|
334
|
|
|
|
|
|
|
335
|
|
|
|
357
|
pool = Pool(processes=self.Num_Chn)
|
|
336
|
pool = Pool(processes=self.Num_Chn)
|
|
|
358
|
args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
|
|
337
|
args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
|
|
|
359
|
objs = [self for __ in range(self.Num_Chn)]
|
|
338
|
objs = [self for __ in range(self.Num_Chn)]
|
|
|
360
|
attrs = list(zip(objs, args))
|
|
339
|
attrs = list(zip(objs, args))
|
|
|
361
|
DGauFitParam = pool.map(target, attrs)
|
|
340
|
gauSPC = pool.map(target, attrs)
|
|
|
362
|
# Parameters:
|
|
341
|
dataOut.SPCparam = numpy.asarray(SPCparam)
|
|
|
363
|
# 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
|
|
342
|
|
|
|
364
|
dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
|
|
343
|
''' Parameters:
|
|
|
365
|
|
|
344
|
1. Amplitude
|
|
|
366
|
# Double Gaussian Curves
|
|
345
|
2. Shift
|
|
|
367
|
gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
|
|
346
|
3. Width
|
|
|
368
|
gau0[:] = numpy.NaN
|
|
347
|
4. Power
|
|
|
369
|
gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
|
|
348
|
'''
|
|
|
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
|
|
|
|
|
|
394
|
|
|
349
|
|
|
|
395
|
def FitGau(self, X):
|
|
350
|
def FitGau(self, X):
|
|
|
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
|
|
|
351
|
|
|
|
407
|
# print ('stop 0')
|
|
352
|
Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
|
|
|
|
|
|
353
|
|
|
|
|
|
|
354
|
SPCparam = []
|
|
|
|
|
|
355
|
SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
|
|
|
|
|
|
356
|
SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
|
|
|
|
|
|
357
|
SPC_ch1[:] = 0#numpy.NaN
|
|
|
|
|
|
358
|
SPC_ch2[:] = 0#numpy.NaN
|
|
|
|
|
|
359
|
|
|
|
408
|
|
|
360
|
|
|
|
409
|
# 5 parameters, 2 Gaussians
|
|
|
|
|
|
410
|
DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
|
|
|
|
|
|
411
|
DGauFitParam[:] = numpy.NaN
|
|
|
|
|
|
412
|
|
|
361
|
|
|
|
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')
|
|
|
|
|
|
419
|
for ht in range(self.Num_Hei):
|
|
362
|
for ht in range(self.Num_Hei):
|
|
|
420
|
# print (ht)
|
|
363
|
|
|
|
421
|
# print ('stop 2')
|
|
364
|
|
|
|
422
|
# Spectra at each range
|
|
|
|
|
|
423
|
spc = numpy.asarray(self.spc)[ch,:,ht]
|
|
365
|
spc = numpy.asarray(self.spc)[ch,:,ht]
|
|
|
424
|
snr = ( spc.mean() - wnoise ) / wnoise
|
|
|
|
|
|
425
|
snrdB = 10.*numpy.log10(snr)
|
|
|
|
|
|
426
|
|
|
366
|
|
|
|
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')
|
|
|
|
|
|
437
|
#############################################
|
|
367
|
#############################################
|
|
|
438
|
# normalizing spc and noise
|
|
368
|
# normalizing spc and noise
|
|
|
439
|
# This part differs from gg1
|
|
369
|
# This part differs from gg1
|
|
|
440
|
# spc_norm_max = max(spc) #commented by D. Scipión 19.03.2021
|
|
370
|
spc_norm_max = max(spc)
|
|
|
441
|
#spc = spc / spc_norm_max
|
|
371
|
#spc = spc / spc_norm_max
|
|
|
442
|
# pnoise = pnoise #/ spc_norm_max #commented by D. Scipión 19.03.2021
|
|
372
|
pnoise = pnoise #/ spc_norm_max
|
|
|
443
|
#############################################
|
|
373
|
#############################################
|
|
|
444
|
|
|
374
|
|
|
|
445
|
# print ('stop 2.1')
|
|
|
|
|
|
446
|
fatspectra=1.0
|
|
375
|
fatspectra=1.0
|
|
|
447
|
# noise per channel.... we might want to use the noise at each range
|
|
|
|
|
|
448
|
|
|
376
|
|
|
|
449
|
# wnoise = noise_ #/ spc_norm_max #commented by D. Scipión 19.03.2021
|
|
377
|
wnoise = noise_ #/ spc_norm_max
|
|
|
450
|
#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
|
|
378
|
#wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
|
|
|
451
|
#if wnoise>1.1*pnoise: # to be tested later
|
|
379
|
#if wnoise>1.1*pnoise: # to be tested later
|
|
|
452
|
# wnoise=pnoise
|
|
380
|
# wnoise=pnoise
|
|
|
453
|
# noisebl = wnoise*0.9
|
|
381
|
noisebl=wnoise*0.9;
|
|
|
454
|
# noisebh = wnoise*1.1
|
|
382
|
noisebh=wnoise*1.1
|
|
|
455
|
spc = spc - wnoise # signal
|
|
383
|
spc=spc-wnoise
|
|
|
456
|
|
|
384
|
|
|
|
457
|
# print ('stop 2.2')
|
|
|
|
|
|
458
|
minx=numpy.argmin(spc)
|
|
385
|
minx=numpy.argmin(spc)
|
|
|
459
|
#spcs=spc.copy()
|
|
386
|
#spcs=spc.copy()
|
|
|
460
|
spcs=numpy.roll(spc,-minx)
|
|
387
|
spcs=numpy.roll(spc,-minx)
|
|
|
461
|
cum=numpy.cumsum(spcs)
|
|
388
|
cum=numpy.cumsum(spcs)
|
|
|
462
|
# tot_noise = wnoise * self.Num_Bin #64;
|
|
389
|
tot_noise=wnoise * self.Num_Bin #64;
|
|
|
463
|
|
|
390
|
|
|
|
464
|
# print ('stop 2.3')
|
|
391
|
snr = sum(spcs)/tot_noise
|
|
|
465
|
# snr = sum(spcs) / tot_noise
|
|
392
|
snrdB=10.*numpy.log10(snr)
|
|
|
466
|
# snrdB = 10.*numpy.log10(snr)
|
|
393
|
|
|
|
467
|
#print ('stop 3')
|
|
394
|
if snrdB < SNRlimit :
|
|
|
468
|
# if snrdB < SNRlimit :
|
|
395
|
snr = numpy.NaN
|
|
|
469
|
# snr = numpy.NaN
|
|
396
|
SPC_ch1[:,ht] = 0#numpy.NaN
|
|
|
470
|
# SPC_ch1[:,ht] = 0#numpy.NaN
|
|
397
|
SPC_ch1[:,ht] = 0#numpy.NaN
|
|
|
471
|
# SPC_ch1[:,ht] = 0#numpy.NaN
|
|
398
|
SPCparam = (SPC_ch1,SPC_ch2)
|
|
|
472
|
# SPCparam = (SPC_ch1,SPC_ch2)
|
|
399
|
continue
|
|
|
473
|
# print ('SNR less than SNRth')
|
|
|
|
|
|
474
|
# continue
|
|
|
|
|
|
475
|
|
|
400
|
|
|
|
476
|
|
|
401
|
|
|
|
477
|
#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
|
|
402
|
#if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
|
|
|
478
|
# return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
|
|
403
|
# return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
|
|
|
479
|
# print ('stop 4')
|
|
404
|
|
|
|
480
|
cummax = max(cum)
|
|
405
|
cummax=max(cum);
|
|
|
481
|
epsi=0.08*fatspectra # cumsum to narrow down the energy region
|
|
406
|
epsi=0.08*fatspectra # cumsum to narrow down the energy region
|
|
|
482
|
cumlo = cummax * epsi
|
|
407
|
cumlo=cummax*epsi;
|
|
|
483
|
cumhi=cummax*(1-epsi)
|
|
408
|
cumhi=cummax*(1-epsi)
|
|
|
484
|
powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
|
|
409
|
powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
|
|
|
485
|
|
|
410
|
|
|
|
486
|
# print ('stop 5')
|
|
411
|
|
|
|
487
|
if len(powerindex) < 1:# case for powerindex 0
|
|
412
|
if len(powerindex) < 1:# case for powerindex 0
|
|
|
488
|
# print ('powerindex < 1')
|
|
|
|
|
|
489
|
continue
|
|
413
|
continue
|
|
|
490
|
powerlo=powerindex[0]
|
|
414
|
powerlo=powerindex[0]
|
|
|
491
|
powerhi=powerindex[-1]
|
|
415
|
powerhi=powerindex[-1]
|
|
|
492
|
powerwidth=powerhi-powerlo
|
|
416
|
powerwidth=powerhi-powerlo
|
|
|
493
|
if powerwidth <= 1:
|
|
|
|
|
|
494
|
# print('powerwidth <= 1')
|
|
|
|
|
|
495
|
continue
|
|
|
|
|
|
496
|
|
|
417
|
|
|
|
497
|
# print ('stop 6')
|
|
|
|
|
|
498
|
firstpeak=powerlo+powerwidth/10.# first gaussian energy location
|
|
418
|
firstpeak=powerlo+powerwidth/10.# first gaussian energy location
|
|
|
499
|
secondpeak=powerhi-powerwidth/10.#second gaussian energy location
|
|
419
|
secondpeak=powerhi-powerwidth/10.#second gaussian energy location
|
|
|
500
|
midpeak=(firstpeak+secondpeak)/2.
|
|
420
|
midpeak=(firstpeak+secondpeak)/2.
|
|
|
@@
-502,6
+422,7
class GaussianFit(Operation):
|
|
|
502
|
secondamp=spcs[int(secondpeak)]
|
|
422
|
secondamp=spcs[int(secondpeak)]
|
|
|
503
|
midamp=spcs[int(midpeak)]
|
|
423
|
midamp=spcs[int(midpeak)]
|
|
|
504
|
|
|
424
|
|
|
|
|
|
|
425
|
x=numpy.arange( self.Num_Bin )
|
|
|
505
|
y_data=spc+wnoise
|
|
426
|
y_data=spc+wnoise
|
|
|
506
|
|
|
427
|
|
|
|
507
|
''' single Gaussian '''
|
|
428
|
''' single Gaussian '''
|
|
|
@@
-510,14
+431,12
class GaussianFit(Operation):
|
|
|
510
|
power0=2.
|
|
431
|
power0=2.
|
|
|
511
|
amplitude0=midamp
|
|
432
|
amplitude0=midamp
|
|
|
512
|
state0=[shift0,width0,amplitude0,power0,wnoise]
|
|
433
|
state0=[shift0,width0,amplitude0,power0,wnoise]
|
|
|
513
|
bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
|
|
434
|
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)
|
|
435
|
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
|
|
|
436
|
|
|
|
518
|
chiSq1=lsq1[1]
|
|
437
|
chiSq1=lsq1[1];
|
|
|
|
|
|
438
|
|
|
|
519
|
|
|
439
|
|
|
|
520
|
# print ('stop 8')
|
|
|
|
|
|
521
|
if fatspectra<1.0 and powerwidth<4:
|
|
440
|
if fatspectra<1.0 and powerwidth<4:
|
|
|
522
|
choice=0
|
|
441
|
choice=0
|
|
|
523
|
Amplitude0=lsq1[0][2]
|
|
442
|
Amplitude0=lsq1[0][2]
|
|
|
@@
-532,33
+451,30
class GaussianFit(Operation):
|
|
|
532
|
#return (numpy.array([shift0,width0,Amplitude0,p0]),
|
|
451
|
#return (numpy.array([shift0,width0,Amplitude0,p0]),
|
|
|
533
|
# numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
|
|
452
|
# numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
|
|
|
534
|
|
|
453
|
|
|
|
535
|
# print ('stop 9')
|
|
454
|
''' two gaussians '''
|
|
|
536
|
''' two Gaussians '''
|
|
|
|
|
|
537
|
#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
|
|
455
|
#shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
|
|
|
538
|
shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
|
|
456
|
shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
|
|
|
539
|
shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
|
|
457
|
shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
|
|
|
540
|
width0 = powerwidth/6.
|
|
458
|
width0=powerwidth/6.;
|
|
|
541
|
width1=width0
|
|
459
|
width1=width0
|
|
|
542
|
power0 = 2.
|
|
460
|
power0=2.;
|
|
|
543
|
power1=power0
|
|
461
|
power1=power0
|
|
|
544
|
amplitude0 = firstamp
|
|
462
|
amplitude0=firstamp;
|
|
|
545
|
amplitude1=secondamp
|
|
463
|
amplitude1=secondamp
|
|
|
546
|
state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
|
|
464
|
state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
|
|
|
547
|
#bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
|
|
465
|
#bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(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))
|
|
466
|
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))
|
|
|
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))
|
|
467
|
#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))
|
|
|
550
|
|
|
468
|
|
|
|
551
|
# print ('stop 10')
|
|
|
|
|
|
552
|
lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
|
|
469
|
lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
|
|
|
553
|
|
|
470
|
|
|
|
554
|
# print ('stop 11')
|
|
|
|
|
|
555
|
chiSq2 = lsq2[1]
|
|
|
|
|
|
556
|
|
|
471
|
|
|
|
557
|
# print ('stop 12')
|
|
472
|
chiSq2=lsq2[1];
|
|
|
|
|
|
473
|
|
|
|
|
|
|
474
|
|
|
|
558
|
|
|
475
|
|
|
|
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)
|
|
476
|
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)
|
|
|
560
|
|
|
477
|
|
|
|
561
|
# print ('stop 13')
|
|
|
|
|
|
562
|
if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
|
|
478
|
if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
|
|
|
563
|
if oneG:
|
|
479
|
if oneG:
|
|
|
564
|
choice=0
|
|
480
|
choice=0
|
|
|
@@
-566,8
+482,8
class GaussianFit(Operation):
|
|
|
566
|
w1=lsq2[0][1]; w2=lsq2[0][5]
|
|
482
|
w1=lsq2[0][1]; w2=lsq2[0][5]
|
|
|
567
|
a1=lsq2[0][2]; a2=lsq2[0][6]
|
|
483
|
a1=lsq2[0][2]; a2=lsq2[0][6]
|
|
|
568
|
p1=lsq2[0][3]; p2=lsq2[0][7]
|
|
484
|
p1=lsq2[0][3]; p2=lsq2[0][7]
|
|
|
569
|
s1 = (2**(1+1./p1))*scipy.special.gamma(1./p1)/p1
|
|
485
|
s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
|
|
|
570
|
s2 = (2**(1+1./p2))*scipy.special.gamma(1./p2)/p2
|
|
486
|
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
|
|
487
|
gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
|
|
|
572
|
|
|
488
|
|
|
|
573
|
if gp1>gp2:
|
|
489
|
if gp1>gp2:
|
|
|
@@
-588,19
+504,16
class GaussianFit(Operation):
|
|
|
588
|
else: # with low SNR go to the most energetic peak
|
|
504
|
else: # with low SNR go to the most energetic peak
|
|
|
589
|
choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
|
|
505
|
choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
|
|
|
590
|
|
|
506
|
|
|
|
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
|
|
|
507
|
|
|
|
597
|
# max_vel = 1.0
|
|
508
|
shift0=lsq2[0][0];
|
|
|
598
|
# Va = max(Vrange)
|
|
509
|
vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
|
|
|
599
|
# deltav = Vrange[1]-Vrange[0]
|
|
510
|
shift1=lsq2[0][4];
|
|
|
600
|
# print ('stop 15')
|
|
511
|
vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
|
|
|
|
|
|
512
|
|
|
|
|
|
|
513
|
max_vel = 1.0
|
|
|
|
|
|
514
|
|
|
|
601
|
#first peak will be 0, second peak will be 1
|
|
515
|
#first peak will be 0, second peak will be 1
|
|
|
602
|
# if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.Scipión 19.03.2021
|
|
516
|
if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
|
|
|
603
|
if vel0 > -Va and vel0 < Va : #first peak is in the correct range
|
|
|
|
|
|
604
|
shift0=lsq2[0][0]
|
|
517
|
shift0=lsq2[0][0]
|
|
|
605
|
width0=lsq2[0][1]
|
|
518
|
width0=lsq2[0][1]
|
|
|
606
|
Amplitude0=lsq2[0][2]
|
|
519
|
Amplitude0=lsq2[0][2]
|
|
|
@@
-624,47
+537,38
class GaussianFit(Operation):
|
|
|
624
|
noise=lsq2[0][8]
|
|
537
|
noise=lsq2[0][8]
|
|
|
625
|
|
|
538
|
|
|
|
626
|
if Amplitude0<0.05: # in case the peak is noise
|
|
539
|
if Amplitude0<0.05: # in case the peak is noise
|
|
|
627
|
shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
|
|
540
|
shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
|
|
|
628
|
if Amplitude1<0.05:
|
|
541
|
if Amplitude1<0.05:
|
|
|
629
|
shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
|
|
542
|
shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
|
|
|
630
|
|
|
543
|
|
|
|
631
|
# print ('stop 16 ')
|
|
544
|
|
|
|
632
|
# SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
|
|
545
|
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)
|
|
546
|
SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
|
|
|
634
|
# SPCparam = (SPC_ch1,SPC_ch2)
|
|
547
|
SPCparam = (SPC_ch1,SPC_ch2)
|
|
|
635
|
|
|
548
|
|
|
|
636
|
DGauFitParam[0,ht,0] = noise
|
|
549
|
|
|
|
637
|
DGauFitParam[0,ht,1] = noise
|
|
550
|
return GauSPC
|
|
|
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
|
|
|
|
|
|
652
|
|
|
551
|
|
|
|
653
|
def y_model1(self,x,state):
|
|
552
|
def y_model1(self,x,state):
|
|
|
654
|
shift0,width0,amplitude0,power0,noise=state
|
|
553
|
shift0,width0,amplitude0,power0,noise=state
|
|
|
655
|
model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
|
|
554
|
model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
|
|
|
|
|
|
555
|
|
|
|
656
|
model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
|
|
556
|
model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
|
|
|
|
|
|
557
|
|
|
|
657
|
model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
|
|
558
|
model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
|
|
|
658
|
return model0+model0u+model0d+noise
|
|
559
|
return model0+model0u+model0d+noise
|
|
|
659
|
|
|
560
|
|
|
|
660
|
def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
|
|
561
|
def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
|
|
|
661
|
shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
|
|
562
|
shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
|
|
|
662
|
model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
|
|
563
|
model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
|
|
|
|
|
|
564
|
|
|
|
663
|
model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
|
|
565
|
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
|
|
|
566
|
|
|
|
|
|
|
567
|
model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
|
|
|
666
|
model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
|
|
568
|
model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
|
|
|
|
|
|
569
|
|
|
|
667
|
model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
|
|
570
|
model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
|
|
|
|
|
|
571
|
|
|
|
668
|
model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
|
|
572
|
model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
|
|
|
669
|
return model0+model0u+model0d+model1+model1u+model1d+noise
|
|
573
|
return model0+model0u+model0d+model1+model1u+model1d+noise
|
|
|
670
|
|
|
574
|
|
|
|
@@
-697,10
+601,31
class PrecipitationProc(Operation):
|
|
|
697
|
Operation.__init__(self)
|
|
601
|
Operation.__init__(self)
|
|
|
698
|
self.i=0
|
|
602
|
self.i=0
|
|
|
699
|
|
|
603
|
|
|
|
|
|
|
604
|
|
|
|
|
|
|
605
|
def gaus(self,xSamples,Amp,Mu,Sigma):
|
|
|
|
|
|
606
|
return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
|
|
|
|
|
|
607
|
|
|
|
|
|
|
608
|
|
|
|
|
|
|
609
|
|
|
|
|
|
|
610
|
def Moments(self, ySamples, xSamples):
|
|
|
|
|
|
611
|
Pot = numpy.nansum( ySamples ) # Potencia, momento 0
|
|
|
|
|
|
612
|
yNorm = ySamples / Pot
|
|
|
|
|
|
613
|
|
|
|
|
|
|
614
|
Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
|
|
|
|
|
|
615
|
Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
|
|
|
|
|
|
616
|
Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
|
|
|
|
|
|
617
|
|
|
|
|
|
|
618
|
return numpy.array([Pot, Vr, Desv])
|
|
|
|
|
|
619
|
|
|
|
700
|
def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
|
|
620
|
def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
|
|
|
701
|
tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30):
|
|
621
|
tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
|
|
|
|
|
|
622
|
|
|
|
702
|
|
|
623
|
|
|
|
703
|
# print ('Entering PrecepitationProc ... ')
|
|
624
|
Velrange = dataOut.spcparam_range[2]
|
|
|
|
|
|
625
|
FrecRange = dataOut.spcparam_range[0]
|
|
|
|
|
|
626
|
|
|
|
|
|
|
627
|
dV= Velrange[1]-Velrange[0]
|
|
|
|
|
|
628
|
dF= FrecRange[1]-FrecRange[0]
|
|
|
704
|
|
|
629
|
|
|
|
705
|
if radar == "MIRA35C" :
|
|
630
|
if radar == "MIRA35C" :
|
|
|
706
|
|
|
631
|
|
|
|
@@
-712,17
+637,18
class PrecipitationProc(Operation):
|
|
|
712
|
|
|
637
|
|
|
|
713
|
else:
|
|
638
|
else:
|
|
|
714
|
|
|
639
|
|
|
|
715
|
self.spc = dataOut.data_pre[0].copy()
|
|
640
|
self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
|
|
|
|
|
|
641
|
|
|
|
|
|
|
642
|
"""NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
|
|
|
716
|
|
|
643
|
|
|
|
717
|
#NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
|
|
|
|
|
|
718
|
self.spc[:,:,0:7]= numpy.NaN
|
|
644
|
self.spc[:,:,0:7]= numpy.NaN
|
|
|
719
|
|
|
645
|
|
|
|
|
|
|
646
|
"""##########################################"""
|
|
|
|
|
|
647
|
|
|
|
720
|
self.Num_Hei = self.spc.shape[2]
|
|
648
|
self.Num_Hei = self.spc.shape[2]
|
|
|
721
|
self.Num_Bin = self.spc.shape[1]
|
|
649
|
self.Num_Bin = self.spc.shape[1]
|
|
|
722
|
self.Num_Chn = self.spc.shape[0]
|
|
650
|
self.Num_Chn = self.spc.shape[0]
|
|
|
723
|
|
|
651
|
|
|
|
724
|
VelRange = dataOut.spc_range[2]
|
|
|
|
|
|
725
|
|
|
|
|
|
|
726
|
''' Se obtiene la constante del RADAR '''
|
|
652
|
''' Se obtiene la constante del RADAR '''
|
|
|
727
|
|
|
653
|
|
|
|
728
|
self.Pt = Pt
|
|
654
|
self.Pt = Pt
|
|
|
@@
-733,71
+659,102
class PrecipitationProc(Operation):
|
|
|
733
|
self.tauW = tauW
|
|
659
|
self.tauW = tauW
|
|
|
734
|
self.ThetaT = ThetaT
|
|
660
|
self.ThetaT = ThetaT
|
|
|
735
|
self.ThetaR = ThetaR
|
|
661
|
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
|
|
|
|
|
|
739
|
|
|
662
|
|
|
|
740
|
Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
|
|
663
|
Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
|
|
|
741
|
Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
|
|
664
|
Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
|
|
|
742
|
RadarConstant = 10e-26 * Numerator / Denominator #
|
|
665
|
RadarConstant = 10e-26 * Numerator / Denominator #
|
|
|
743
|
ExpConstant = 10**(40/10) #Constante Experimental
|
|
|
|
|
|
744
|
|
|
666
|
|
|
|
745
|
SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
|
|
667
|
''' ============================= '''
|
|
|
746
|
for i in range(self.Num_Chn):
|
|
668
|
|
|
|
747
|
SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
|
|
669
|
self.spc[0] = (self.spc[0]-dataOut.noise[0])
|
|
|
748
|
SignalPower[numpy.where(SignalPower < 0)] = 1e-20
|
|
670
|
self.spc[1] = (self.spc[1]-dataOut.noise[1])
|
|
|
749
|
|
|
671
|
self.spc[2] = (self.spc[2]-dataOut.noise[2])
|
|
|
750
|
SPCmean = numpy.mean(SignalPower, 0)
|
|
672
|
|
|
|
751
|
Pr = SPCmean[:,:]/dataOut.normFactor
|
|
673
|
self.spc[ numpy.where(self.spc < 0)] = 0
|
|
|
752
|
|
|
674
|
|
|
|
753
|
# Declaring auxiliary variables
|
|
675
|
SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
|
|
|
754
|
Range = dataOut.heightList*1000. #Range in m
|
|
676
|
SPCmean[ numpy.where(SPCmean < 0)] = 0
|
|
|
755
|
# replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
|
|
677
|
|
|
|
756
|
rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
|
|
678
|
ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
757
|
zMtrx = rMtrx+Altitude
|
|
679
|
ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
758
|
# replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
|
|
680
|
ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
759
|
VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
|
|
681
|
|
|
|
760
|
|
|
682
|
Pr = SPCmean[:,:]
|
|
|
761
|
# height dependence to air density Foote and Du Toit (1969)
|
|
683
|
|
|
|
762
|
delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
|
|
684
|
VelMeteoro = numpy.mean(SPCmean,axis=0)
|
|
|
763
|
VMtrx = VelMtrx / delv_z #Normalized velocity
|
|
685
|
|
|
|
764
|
VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
|
|
686
|
D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
765
|
# Diameter is related to the fall speed of falling drops
|
|
687
|
SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
766
|
D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
|
|
688
|
N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
767
|
# Only valid for D>= 0.16 mm
|
|
689
|
V_mean = numpy.zeros(self.Num_Hei)
|
|
|
768
|
D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
|
|
690
|
del_V = numpy.zeros(self.Num_Hei)
|
|
|
769
|
|
|
691
|
Z = numpy.zeros(self.Num_Hei)
|
|
|
770
|
#Calculate Radar Reflectivity ETAn
|
|
692
|
Ze = numpy.zeros(self.Num_Hei)
|
|
|
771
|
ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
|
|
693
|
RR = numpy.zeros(self.Num_Hei)
|
|
|
772
|
ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
|
|
694
|
|
|
|
773
|
# Radar Cross Section
|
|
695
|
Range = dataOut.heightList*1000.
|
|
|
774
|
sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
|
|
696
|
|
|
|
775
|
# Drop Size Distribution
|
|
697
|
for R in range(self.Num_Hei):
|
|
|
776
|
DSD = ETAn / sigmaD
|
|
698
|
|
|
|
777
|
# Equivalente Reflectivy
|
|
699
|
h = Range[R] + Altitude #Range from ground to radar pulse altitude
|
|
|
778
|
Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
|
|
700
|
del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
|
|
|
779
|
Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
|
|
701
|
|
|
|
780
|
# RainFall Rate
|
|
702
|
D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
|
|
|
781
|
RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
|
|
703
|
|
|
|
782
|
|
|
704
|
'''NOTA: ETA(n) dn = ETA(f) df
|
|
|
783
|
# Censoring the data
|
|
705
|
|
|
|
784
|
# Removing data with SNRth < 0dB se debe considerar el SNR por canal
|
|
706
|
dn = 1 Diferencial de muestreo
|
|
|
785
|
SNRth = 10**(SNRdBlimit/10) #-30dB
|
|
707
|
df = ETA(n) / ETA(f)
|
|
|
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
|
|
708
|
|
|
|
787
|
W = numpy.nanmean(dataOut.data_dop,0)
|
|
709
|
'''
|
|
|
788
|
W[novalid] = numpy.NaN
|
|
710
|
|
|
|
789
|
Ze_org[novalid] = numpy.NaN
|
|
711
|
ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
|
|
|
790
|
RR[novalid] = numpy.NaN
|
|
712
|
|
|
|
|
|
|
713
|
ETAv[:,R]=ETAn[:,R]/dV
|
|
|
|
|
|
714
|
|
|
|
|
|
|
715
|
ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
|
|
|
|
|
|
716
|
|
|
|
|
|
|
717
|
SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
|
|
|
|
|
|
718
|
|
|
|
|
|
|
719
|
N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
|
|
|
|
|
|
720
|
|
|
|
|
|
|
721
|
DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
|
|
|
|
|
|
722
|
|
|
|
|
|
|
723
|
try:
|
|
|
|
|
|
724
|
popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
|
|
|
|
|
|
725
|
except:
|
|
|
|
|
|
726
|
popt01=numpy.zeros(3)
|
|
|
|
|
|
727
|
popt01[1]= DMoments[1]
|
|
|
|
|
|
728
|
|
|
|
|
|
|
729
|
if popt01[1]<0 or popt01[1]>20:
|
|
|
|
|
|
730
|
popt01[1]=numpy.NaN
|
|
|
|
|
|
731
|
|
|
|
|
|
|
732
|
|
|
|
|
|
|
733
|
V_mean[R]=popt01[1]
|
|
|
|
|
|
734
|
|
|
|
|
|
|
735
|
Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
|
|
|
|
|
|
736
|
|
|
|
|
|
|
737
|
RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
|
|
|
|
|
|
738
|
|
|
|
|
|
|
739
|
Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
|
|
|
|
|
|
740
|
|
|
|
|
|
|
741
|
|
|
|
|
|
|
742
|
|
|
|
|
|
|
743
|
RR2 = (Z/200)**(1/1.6)
|
|
|
|
|
|
744
|
dBRR = 10*numpy.log10(RR)
|
|
|
|
|
|
745
|
dBRR2 = 10*numpy.log10(RR2)
|
|
|
|
|
|
746
|
|
|
|
|
|
|
747
|
dBZe = 10*numpy.log10(Ze)
|
|
|
|
|
|
748
|
dBZ = 10*numpy.log10(Z)
|
|
|
791
|
|
|
749
|
|
|
|
792
|
dataOut.data_output = RR[8]
|
|
750
|
dataOut.data_output = RR[8]
|
|
|
793
|
dataOut.data_param = numpy.ones([3,self.Num_Hei])
|
|
751
|
dataOut.data_param = numpy.ones([3,self.Num_Hei])
|
|
|
794
|
dataOut.channelList = [0,1,2]
|
|
752
|
dataOut.channelList = [0,1,2]
|
|
|
795
|
|
|
753
|
|
|
|
796
|
dataOut.data_param[0]=10*numpy.log10(Ze_org)
|
|
754
|
dataOut.data_param[0]=dBZ
|
|
|
797
|
dataOut.data_param[1]=-W
|
|
755
|
dataOut.data_param[1]=V_mean
|
|
|
798
|
dataOut.data_param[2]=RR
|
|
756
|
dataOut.data_param[2]=RR
|
|
|
799
|
|
|
757
|
|
|
|
800
|
# print ('Leaving PrecepitationProc ... ')
|
|
|
|
|
|
801
|
return dataOut
|
|
758
|
return dataOut
|
|
|
802
|
|
|
759
|
|
|
|
803
|
def dBZeMODE2(self, dataOut): # Processing for MIRA35C
|
|
760
|
def dBZeMODE2(self, dataOut): # Processing for MIRA35C
|
|
|
@@
-814,7
+771,7
class PrecipitationProc(Operation):
|
|
|
814
|
|
|
771
|
|
|
|
815
|
ETA = numpy.sum(SNR,1)
|
|
772
|
ETA = numpy.sum(SNR,1)
|
|
|
816
|
|
|
773
|
|
|
|
817
|
ETA = numpy.where(ETA != 0. , ETA, numpy.NaN)
|
|
774
|
ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
|
|
|
818
|
|
|
775
|
|
|
|
819
|
Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
|
|
776
|
Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
|
|
|
820
|
|
|
777
|
|
|
|
@@
-852,46
+809,40
class PrecipitationProc(Operation):
|
|
|
852
|
class FullSpectralAnalysis(Operation):
|
|
809
|
class FullSpectralAnalysis(Operation):
|
|
|
853
|
|
|
810
|
|
|
|
854
|
"""
|
|
811
|
"""
|
|
|
855
|
Function that implements Full Spectral Analysis technique.
|
|
812
|
Function that implements Full Spectral Analisys technique.
|
|
|
856
|
|
|
813
|
|
|
|
857
|
Input:
|
|
814
|
Input:
|
|
|
858
|
self.dataOut.data_pre : SelfSpectra and CrossSpectra data
|
|
815
|
self.dataOut.data_pre : SelfSpectra and CrossSPectra data
|
|
|
859
|
self.dataOut.groupList : Pairlist of channels
|
|
816
|
self.dataOut.groupList : Pairlist of channels
|
|
|
860
|
self.dataOut.ChanDist : Physical distance between receivers
|
|
817
|
self.dataOut.ChanDist : Physical distance between receivers
|
|
|
861
|
|
|
818
|
|
|
|
862
|
|
|
819
|
|
|
|
863
|
Output:
|
|
820
|
Output:
|
|
|
864
|
|
|
821
|
|
|
|
865
|
self.dataOut.data_output : Zonal wind, Meridional wind, and Vertical wind
|
|
822
|
self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
|
|
|
866
|
|
|
823
|
|
|
|
867
|
|
|
824
|
|
|
|
868
|
Parameters affected: Winds, height range, SNR
|
|
825
|
Parameters affected: Winds, height range, SNR
|
|
|
869
|
|
|
826
|
|
|
|
870
|
"""
|
|
827
|
"""
|
|
|
871
|
def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
|
|
828
|
def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7):
|
|
|
872
|
minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
|
|
829
|
|
|
|
|
|
|
830
|
self.indice=int(numpy.random.rand()*1000)
|
|
|
873
|
|
|
831
|
|
|
|
874
|
spc = dataOut.data_pre[0].copy()
|
|
832
|
spc = dataOut.data_pre[0].copy()
|
|
|
875
|
cspc = dataOut.data_pre[1]
|
|
833
|
cspc = dataOut.data_pre[1]
|
|
|
876
|
nHeights = spc.shape[2]
|
|
|
|
|
|
877
|
|
|
834
|
|
|
|
878
|
# first_height = 0.75 #km (ref: data header 20170822)
|
|
835
|
"""NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
|
|
|
879
|
# resolution_height = 0.075 #km
|
|
836
|
|
|
|
880
|
'''
|
|
837
|
SNRspc = spc.copy()
|
|
|
881
|
finding height range. check this when radar parameters are changed!
|
|
838
|
SNRspc[:,:,0:7]= numpy.NaN
|
|
|
882
|
'''
|
|
839
|
|
|
|
883
|
if maxheight is not None:
|
|
840
|
"""##########################################"""
|
|
|
884
|
# range_max = math.ceil((maxheight - first_height) / resolution_height) # theoretical
|
|
841
|
|
|
|
885
|
range_max = math.ceil(13.26 * maxheight - 3) # empirical, works better
|
|
842
|
|
|
|
886
|
else:
|
|
843
|
nChannel = spc.shape[0]
|
|
|
887
|
range_max = nHeights
|
|
844
|
nProfiles = spc.shape[1]
|
|
|
888
|
if minheight is not None:
|
|
845
|
nHeights = spc.shape[2]
|
|
|
889
|
# range_min = int((minheight - first_height) / resolution_height) # theoretical
|
|
|
|
|
|
890
|
range_min = int(13.26 * minheight - 5) # empirical, works better
|
|
|
|
|
|
891
|
if range_min < 0:
|
|
|
|
|
|
892
|
range_min = 0
|
|
|
|
|
|
893
|
else:
|
|
|
|
|
|
894
|
range_min = 0
|
|
|
|
|
|
895
|
|
|
846
|
|
|
|
896
|
pairsList = dataOut.groupList
|
|
847
|
pairsList = dataOut.groupList
|
|
|
897
|
if dataOut.ChanDist is not None :
|
|
848
|
if dataOut.ChanDist is not None :
|
|
|
@@
-899,301
+850,329
class FullSpectralAnalysis(Operation):
|
|
|
899
|
else:
|
|
850
|
else:
|
|
|
900
|
ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
|
|
851
|
ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
|
|
|
901
|
|
|
852
|
|
|
|
902
|
# 4 variables: zonal, meridional, vertical, and average SNR
|
|
853
|
FrecRange = dataOut.spc_range[0]
|
|
|
903
|
data_param = numpy.zeros([4,nHeights]) * numpy.NaN
|
|
854
|
|
|
|
904
|
velocityX = numpy.zeros([nHeights]) * numpy.NaN
|
|
855
|
ySamples=numpy.ones([nChannel,nProfiles])
|
|
|
905
|
velocityY = numpy.zeros([nHeights]) * numpy.NaN
|
|
856
|
phase=numpy.ones([nChannel,nProfiles])
|
|
|
906
|
velocityZ = numpy.zeros([nHeights]) * numpy.NaN
|
|
857
|
CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
|
|
|
|
|
|
858
|
coherence=numpy.ones([nChannel,nProfiles])
|
|
|
|
|
|
859
|
PhaseSlope=numpy.ones(nChannel)
|
|
|
|
|
|
860
|
PhaseInter=numpy.ones(nChannel)
|
|
|
|
|
|
861
|
data_SNR=numpy.zeros([nProfiles])
|
|
|
|
|
|
862
|
|
|
|
|
|
|
863
|
data = dataOut.data_pre
|
|
|
|
|
|
864
|
noise = dataOut.noise
|
|
|
|
|
|
865
|
|
|
|
|
|
|
866
|
dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
|
|
|
|
|
|
867
|
|
|
|
|
|
|
868
|
dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
|
|
|
907
|
|
|
869
|
|
|
|
908
|
dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
|
|
|
|
|
|
909
|
|
|
870
|
|
|
|
910
|
'''***********************************************WIND ESTIMATION**************************************'''
|
|
871
|
data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
|
|
|
|
|
|
872
|
|
|
|
|
|
|
873
|
velocityX=[]
|
|
|
|
|
|
874
|
velocityY=[]
|
|
|
|
|
|
875
|
velocityV=[]
|
|
|
|
|
|
876
|
PhaseLine=[]
|
|
|
|
|
|
877
|
|
|
|
|
|
|
878
|
dbSNR = 10*numpy.log10(dataOut.data_SNR)
|
|
|
|
|
|
879
|
dbSNR = numpy.average(dbSNR,0)
|
|
|
|
|
|
880
|
|
|
|
911
|
for Height in range(nHeights):
|
|
881
|
for Height in range(nHeights):
|
|
|
912
|
|
|
882
|
|
|
|
913
|
if Height >= range_min and Height < range_max:
|
|
883
|
[Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
|
|
|
914
|
# error_code will be useful in future analysis
|
|
884
|
PhaseLine = numpy.append(PhaseLine, PhaseSlope)
|
|
|
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
|
|
|
885
|
|
|
|
918
|
if abs(Vzon) < 100. and abs(Vmer) < 100.:
|
|
886
|
if abs(Vzon)<100. and abs(Vzon)> 0.:
|
|
|
919
|
velocityX[Height] = Vzon
|
|
887
|
velocityX=numpy.append(velocityX, Vzon)#Vmag
|
|
|
920
|
velocityY[Height] = -Vmer
|
|
|
|
|
|
921
|
velocityZ[Height] = Vver
|
|
|
|
|
|
922
|
|
|
888
|
|
|
|
923
|
# Censoring data with SNR threshold
|
|
889
|
else:
|
|
|
924
|
dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
|
|
890
|
velocityX=numpy.append(velocityX, numpy.NaN)
|
|
|
|
|
|
891
|
|
|
|
|
|
|
892
|
if abs(Vmer)<100. and abs(Vmer) > 0.:
|
|
|
|
|
|
893
|
velocityY=numpy.append(velocityY, -Vmer)#Vang
|
|
|
|
|
|
894
|
|
|
|
|
|
|
895
|
else:
|
|
|
|
|
|
896
|
velocityY=numpy.append(velocityY, numpy.NaN)
|
|
|
|
|
|
897
|
|
|
|
|
|
|
898
|
if dbSNR[Height] > SNRlimit:
|
|
|
|
|
|
899
|
velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
|
|
|
|
|
|
900
|
else:
|
|
|
|
|
|
901
|
velocityV=numpy.append(velocityV, numpy.NaN)
|
|
|
|
|
|
902
|
|
|
|
|
|
|
903
|
|
|
|
|
|
|
904
|
|
|
|
|
|
|
905
|
'''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR'''
|
|
|
|
|
|
906
|
data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
|
|
|
|
|
|
907
|
data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
|
|
|
|
|
|
908
|
data_output[2] = velocityV#FirstMoment
|
|
|
|
|
|
909
|
|
|
|
|
|
|
910
|
xFrec=FrecRange[0:spc.shape[1]]
|
|
|
|
|
|
911
|
|
|
|
|
|
|
912
|
dataOut.data_output=data_output
|
|
|
925
|
|
|
913
|
|
|
|
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
|
|
|
|
|
|
931
|
return dataOut
|
|
914
|
return dataOut
|
|
|
932
|
|
|
915
|
|
|
|
|
|
|
916
|
|
|
|
933
|
def moving_average(self,x, N=2):
|
|
917
|
def moving_average(self,x, N=2):
|
|
|
934
|
""" convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
|
|
|
|
|
|
935
|
return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
|
|
918
|
return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
|
|
|
936
|
|
|
919
|
|
|
|
937
|
def gaus(self,xSamples,Amp,Mu,Sigma):
|
|
920
|
def gaus(self,xSamples,Amp,Mu,Sigma):
|
|
|
938
|
return Amp * numpy.exp(-0.5*((xSamples - Mu)/Sigma)**2)
|
|
921
|
return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
|
|
|
939
|
|
|
922
|
|
|
|
940
|
def Moments(self, ySamples, xSamples):
|
|
|
|
|
|
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])
|
|
|
|
|
|
947
|
|
|
|
|
|
|
948
|
def StopWindEstimation(self, error_code):
|
|
|
|
|
|
949
|
Vzon = numpy.NaN
|
|
|
|
|
|
950
|
Vmer = numpy.NaN
|
|
|
|
|
|
951
|
Vver = numpy.NaN
|
|
|
|
|
|
952
|
return Vzon, Vmer, Vver, error_code
|
|
|
|
|
|
953
|
|
|
923
|
|
|
|
954
|
def AntiAliasing(self, interval, maxstep):
|
|
|
|
|
|
955
|
"""
|
|
|
|
|
|
956
|
function to prevent errors from aliased values when computing phaseslope
|
|
|
|
|
|
957
|
"""
|
|
|
|
|
|
958
|
antialiased = numpy.zeros(len(interval))
|
|
|
|
|
|
959
|
copyinterval = interval.copy()
|
|
|
|
|
|
960
|
|
|
|
|
|
|
961
|
antialiased[0] = copyinterval[0]
|
|
|
|
|
|
962
|
|
|
|
|
|
|
963
|
for i in range(1,len(antialiased)):
|
|
|
|
|
|
964
|
step = interval[i] - interval[i-1]
|
|
|
|
|
|
965
|
if step > maxstep:
|
|
|
|
|
|
966
|
copyinterval -= 2*numpy.pi
|
|
|
|
|
|
967
|
antialiased[i] = copyinterval[i]
|
|
|
|
|
|
968
|
elif step < maxstep*(-1):
|
|
|
|
|
|
969
|
copyinterval += 2*numpy.pi
|
|
|
|
|
|
970
|
antialiased[i] = copyinterval[i]
|
|
|
|
|
|
971
|
else:
|
|
|
|
|
|
972
|
antialiased[i] = copyinterval[i].copy()
|
|
|
|
|
|
973
|
|
|
924
|
|
|
|
974
|
return antialiased
|
|
925
|
def Moments(self, ySamples, xSamples):
|
|
|
|
|
|
926
|
Pot = numpy.nansum( ySamples ) # Potencia, momento 0
|
|
|
|
|
|
927
|
yNorm = ySamples / Pot
|
|
|
|
|
|
928
|
Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
|
|
|
|
|
|
929
|
Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
|
|
|
|
|
|
930
|
Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
|
|
|
975
|
|
|
931
|
|
|
|
976
|
def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit, NegativeLimit, PositiveLimit, radfreq):
|
|
932
|
return numpy.array([Pot, Vr, Desv])
|
|
|
977
|
"""
|
|
|
|
|
|
978
|
Function that Calculates Zonal, Meridional and Vertical wind velocities.
|
|
|
|
|
|
979
|
Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
|
|
|
|
|
|
980
|
|
|
933
|
|
|
|
981
|
Input:
|
|
934
|
def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
|
|
|
982
|
spc, cspc : self spectra and cross spectra data. In Briggs notation something like S_i*(S_i)_conj, (S_j)_conj respectively.
|
|
|
|
|
|
983
|
pairsList : Pairlist of channels
|
|
|
|
|
|
984
|
ChanDist : array of xi_ij and eta_ij
|
|
|
|
|
|
985
|
Height : height at which data is processed
|
|
|
|
|
|
986
|
noise : noise in [channels] format for specific height
|
|
|
|
|
|
987
|
Abbsisarange : range of the frequencies or velocities
|
|
|
|
|
|
988
|
dbSNR, SNRlimit : signal to noise ratio in db, lower limit
|
|
|
|
|
|
989
|
|
|
935
|
|
|
|
990
|
Output:
|
|
|
|
|
|
991
|
Vzon, Vmer, Vver : wind velocities
|
|
|
|
|
|
992
|
error_code : int that states where code is terminated
|
|
|
|
|
|
993
|
|
|
|
|
|
|
994
|
0 : no error detected
|
|
|
|
|
|
995
|
1 : Gaussian of mean spc exceeds widthlimit
|
|
|
|
|
|
996
|
2 : no Gaussian of mean spc found
|
|
|
|
|
|
997
|
3 : SNR to low or velocity to high -> prec. e.g.
|
|
|
|
|
|
998
|
4 : at least one Gaussian of cspc exceeds widthlimit
|
|
|
|
|
|
999
|
5 : zero out of three cspc Gaussian fits converged
|
|
|
|
|
|
1000
|
6 : phase slope fit could not be found
|
|
|
|
|
|
1001
|
7 : arrays used to fit phase have different length
|
|
|
|
|
|
1002
|
8 : frequency range is either too short (len <= 5) or very long (> 30% of cspc)
|
|
|
|
|
|
1003
|
|
|
936
|
|
|
|
1004
|
"""
|
|
|
|
|
|
1005
|
|
|
937
|
|
|
|
1006
|
error_code = 0
|
|
938
|
ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
|
1007
|
|
|
939
|
phase=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
|
1008
|
nChan = spc.shape[0]
|
|
940
|
CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
|
|
|
1009
|
nProf = spc.shape[1]
|
|
941
|
coherence=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
|
1010
|
nPair = cspc.shape[0]
|
|
942
|
PhaseSlope=numpy.zeros(spc.shape[0])
|
|
|
1011
|
|
|
943
|
PhaseInter=numpy.ones(spc.shape[0])
|
|
|
1012
|
SPC_Samples = numpy.zeros([nChan, nProf]) # for normalized spc values for one height
|
|
944
|
xFrec=AbbsisaRange[0][0:spc.shape[1]]
|
|
|
1013
|
CSPC_Samples = numpy.zeros([nPair, nProf], dtype=numpy.complex_) # for normalized cspc values
|
|
945
|
xVel =AbbsisaRange[2][0:spc.shape[1]]
|
|
|
1014
|
phase = numpy.zeros([nPair, nProf]) # phase between channels
|
|
946
|
Vv=numpy.empty(spc.shape[2])*0
|
|
|
1015
|
PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
|
|
947
|
SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
|
|
|
1016
|
PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
|
|
948
|
|
|
|
1017
|
xFrec = AbbsisaRange[0][:-1] # frequency range
|
|
949
|
SPCmoments = self.Moments(SPCav[:,Height], xVel )
|
|
|
1018
|
xVel = AbbsisaRange[2][:-1] # velocity range
|
|
950
|
CSPCmoments = []
|
|
|
1019
|
xSamples = xFrec # the frequency range is taken
|
|
951
|
cspcNoise = numpy.empty(3)
|
|
|
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]
|
|
|
|
|
|
1029
|
|
|
952
|
|
|
|
1030
|
'''Getting Eij and Nij'''
|
|
953
|
'''Getting Eij and Nij'''
|
|
|
1031
|
Xi01, Xi02, Xi12 = ChanDist[:,0]
|
|
|
|
|
|
1032
|
Eta01, Eta02, Eta12 = ChanDist[:,1]
|
|
|
|
|
|
1033
|
|
|
|
|
|
|
1034
|
# spwd limit - updated by D. Scipión 30.03.2021
|
|
|
|
|
|
1035
|
widthlimit = 10
|
|
|
|
|
|
1036
|
'''************************* SPC is normalized ********************************'''
|
|
|
|
|
|
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)
|
|
|
|
|
|
1042
|
|
|
954
|
|
|
|
1043
|
'''********************** FITTING MEAN SPC GAUSSIAN **********************'''
|
|
955
|
Xi01=ChanDist[0][0]
|
|
|
|
|
|
956
|
Eta01=ChanDist[0][1]
|
|
|
1044
|
|
|
957
|
|
|
|
1045
|
""" the gaussian of the mean: first subtract noise, then normalize. this is legal because
|
|
958
|
Xi02=ChanDist[1][0]
|
|
|
1046
|
you only fit the curve and don't need the absolute value of height for calculation,
|
|
959
|
Eta02=ChanDist[1][1]
|
|
|
1047
|
only for estimation of width. for normalization of cross spectra, you need initial,
|
|
|
|
|
|
1048
|
unnormalized self-spectra With noise.
|
|
|
|
|
|
1049
|
|
|
960
|
|
|
|
1050
|
Technically, you don't even need to normalize the self-spectra, as you only need the
|
|
961
|
Xi12=ChanDist[2][0]
|
|
|
1051
|
width of the peak. However, it was left this way. Note that the normalization has a flaw:
|
|
962
|
Eta12=ChanDist[2][1]
|
|
|
1052
|
due to subtraction of the noise, some values are below zero. Raw "spc" values should be
|
|
|
|
|
|
1053
|
>= 0, as it is the modulus squared of the signals (complex * it's conjugate)
|
|
|
|
|
|
1054
|
"""
|
|
|
|
|
|
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
|
|
|
963
|
|
|
|
1062
|
# Gauss Fit SPC in frequency domain
|
|
964
|
z = spc.copy()
|
|
|
1063
|
if dbSNR > SNRlimit: # only if SNR > SNRth
|
|
965
|
z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
|
|
|
1064
|
try:
|
|
|
|
|
|
1065
|
popt,pcov = curve_fit(self.gaus,xSamples_zoom,SPCMean[xvalid],p0=SPCMoments)
|
|
|
|
|
|
1066
|
if popt[2] <= 0 or popt[2] > widthlimit: # CONDITION
|
|
|
|
|
|
1067
|
return self.StopWindEstimation(error_code = 1)
|
|
|
|
|
|
1068
|
FitGauss = self.gaus(xSamples_zoom,*popt)
|
|
|
|
|
|
1069
|
except :#RuntimeError:
|
|
|
|
|
|
1070
|
return self.StopWindEstimation(error_code = 2)
|
|
|
|
|
|
1071
|
else:
|
|
|
|
|
|
1072
|
return self.StopWindEstimation(error_code = 3)
|
|
|
|
|
|
1073
|
|
|
966
|
|
|
|
1074
|
'''***************************** CSPC Normalization *************************
|
|
967
|
for i in range(spc.shape[0]):
|
|
|
1075
|
The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
|
|
|
|
|
|
1076
|
influence the norm which is not desired. First, a range is identified where the
|
|
|
|
|
|
1077
|
wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
|
|
|
|
|
|
1078
|
around it gets cut off and values replaced by mean determined by the boundary
|
|
|
|
|
|
1079
|
data -> sum_noise (spc is not normalized here, thats why the noise is important)
|
|
|
|
|
|
1080
|
|
|
968
|
|
|
|
1081
|
The sums are then added and multiplied by range/datapoints, because you need
|
|
969
|
'''****** Line of Data SPC ******'''
|
|
|
1082
|
an integral and not a sum for normalization.
|
|
970
|
zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
|
|
|
1083
|
|
|
971
|
|
|
|
1084
|
A norm is found according to Briggs 92.
|
|
972
|
'''****** SPC is normalized ******'''
|
|
|
1085
|
'''
|
|
973
|
SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
|
|
|
1086
|
# for each pair
|
|
974
|
FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
|
|
|
1087
|
for i in range(nPair):
|
|
975
|
|
|
|
1088
|
cspc_norm = cspc[i,:].copy()
|
|
976
|
xSamples = xFrec # Se toma el rango de frecuncias
|
|
|
|
|
|
977
|
ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
|
|
|
|
|
|
978
|
|
|
|
|
|
|
979
|
for i in range(spc.shape[0]):
|
|
|
|
|
|
980
|
|
|
|
|
|
|
981
|
'''****** Line of Data CSPC ******'''
|
|
|
|
|
|
982
|
cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
|
|
|
|
|
|
983
|
SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
|
|
|
|
|
|
984
|
cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
|
|
|
|
|
|
985
|
|
|
|
|
|
|
986
|
'''****** CSPC is normalized with respect to Briggs and Vincent ******'''
|
|
|
1089
|
chan_index0 = pairsList[i][0]
|
|
987
|
chan_index0 = pairsList[i][0]
|
|
|
1090
|
chan_index1 = pairsList[i][1]
|
|
988
|
chan_index1 = pairsList[i][1]
|
|
|
1091
|
CSPC_Samples[i] = cspc_norm / (numpy.sqrt(numpy.nansum(spc_norm[chan_index0])*numpy.nansum(spc_norm[chan_index1])) * delta_x)
|
|
|
|
|
|
1092
|
phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
|
|
|
|
|
|
1093
|
|
|
989
|
|
|
|
1094
|
CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0,xvalid]), xSamples_zoom),
|
|
990
|
CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
|
|
|
1095
|
self.Moments(numpy.abs(CSPC_Samples[1,xvalid]), xSamples_zoom),
|
|
991
|
CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
|
|
|
1096
|
self.Moments(numpy.abs(CSPC_Samples[2,xvalid]), xSamples_zoom)])
|
|
992
|
|
|
|
|
|
|
993
|
CSPCSamples[i] = CSPCNorm
|
|
|
|
|
|
994
|
|
|
|
|
|
|
995
|
coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
|
|
|
|
|
|
996
|
|
|
|
|
|
|
997
|
#coherence[i]= self.moving_average(coherence[i],N=1)
|
|
|
|
|
|
998
|
|
|
|
|
|
|
999
|
phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
|
|
|
|
|
|
1000
|
|
|
|
|
|
|
1001
|
CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
|
|
|
|
|
|
1002
|
self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
|
|
|
|
|
|
1003
|
self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
|
|
|
|
|
|
1004
|
|
|
|
|
|
|
1005
|
|
|
|
|
|
|
1006
|
popt=[1e-10,0,1e-10]
|
|
|
|
|
|
1007
|
popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
|
|
|
|
|
|
1008
|
FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
|
|
|
|
|
|
1009
|
|
|
|
|
|
|
1010
|
CSPCMask01 = numpy.abs(CSPCSamples[0])
|
|
|
|
|
|
1011
|
CSPCMask02 = numpy.abs(CSPCSamples[1])
|
|
|
|
|
|
1012
|
CSPCMask12 = numpy.abs(CSPCSamples[2])
|
|
|
|
|
|
1013
|
|
|
|
|
|
|
1014
|
mask01 = ~numpy.isnan(CSPCMask01)
|
|
|
|
|
|
1015
|
mask02 = ~numpy.isnan(CSPCMask02)
|
|
|
|
|
|
1016
|
mask12 = ~numpy.isnan(CSPCMask12)
|
|
|
|
|
|
1017
|
|
|
|
|
|
|
1018
|
#mask = ~numpy.isnan(CSPCMask01)
|
|
|
|
|
|
1019
|
CSPCMask01 = CSPCMask01[mask01]
|
|
|
|
|
|
1020
|
CSPCMask02 = CSPCMask02[mask02]
|
|
|
|
|
|
1021
|
CSPCMask12 = CSPCMask12[mask12]
|
|
|
|
|
|
1022
|
#CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
|
|
|
1097
|
|
|
1023
|
|
|
|
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))
|
|
|
|
|
|
1100
|
|
|
1024
|
|
|
|
1101
|
'''*******************************FIT GAUSS CSPC************************************'''
|
|
1025
|
|
|
|
|
|
|
1026
|
'''***Fit Gauss CSPC01***'''
|
|
|
|
|
|
1027
|
if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
|
|
|
1102
|
try:
|
|
1028
|
try:
|
|
|
1103
|
popt01,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[0][xvalid]),p0=CSPCmoments[0])
|
|
1029
|
popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
|
|
|
1104
|
if popt01[2] > widthlimit: # CONDITION
|
|
1030
|
popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
|
|
|
1105
|
return self.StopWindEstimation(error_code = 4)
|
|
1031
|
popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
|
|
|
1106
|
popt02,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[1][xvalid]),p0=CSPCmoments[1])
|
|
1032
|
FitGauss01 = self.gaus(xSamples,*popt01)
|
|
|
1107
|
if popt02[2] > widthlimit: # CONDITION
|
|
1033
|
FitGauss02 = self.gaus(xSamples,*popt02)
|
|
|
1108
|
return self.StopWindEstimation(error_code = 4)
|
|
1034
|
FitGauss12 = self.gaus(xSamples,*popt12)
|
|
|
1109
|
popt12,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[2][xvalid]),p0=CSPCmoments[2])
|
|
|
|
|
|
1110
|
if popt12[2] > widthlimit: # CONDITION
|
|
|
|
|
|
1111
|
return self.StopWindEstimation(error_code = 4)
|
|
|
|
|
|
1112
|
|
|
|
|
|
|
1113
|
FitGauss01 = self.gaus(xSamples_zoom, *popt01)
|
|
|
|
|
|
1114
|
FitGauss02 = self.gaus(xSamples_zoom, *popt02)
|
|
|
|
|
|
1115
|
FitGauss12 = self.gaus(xSamples_zoom, *popt12)
|
|
|
|
|
|
1116
|
except:
|
|
1035
|
except:
|
|
|
1117
|
return self.StopWindEstimation(error_code = 5)
|
|
1036
|
FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
|
|
|
|
|
|
1037
|
FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
|
|
|
|
|
|
1038
|
FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
|
|
|
|
|
|
1039
|
|
|
|
|
|
|
1040
|
|
|
|
|
|
|
1041
|
CSPCopt = numpy.vstack([popt01,popt02,popt12])
|
|
|
|
|
|
1042
|
|
|
|
|
|
|
1043
|
'''****** Getting fij width ******'''
|
|
|
|
|
|
1044
|
|
|
|
|
|
|
1045
|
yMean = numpy.average(ySamples, axis=0) # ySamples[0]
|
|
|
|
|
|
1046
|
|
|
|
|
|
|
1047
|
'''******* Getting fitting Gaussian *******'''
|
|
|
|
|
|
1048
|
meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
|
|
|
|
|
|
1049
|
sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
|
|
|
1118
|
|
|
1050
|
|
|
|
|
|
|
1051
|
yMoments = self.Moments(yMean, xSamples)
|
|
|
1119
|
|
|
1052
|
|
|
|
1120
|
'''************* Getting Fij ***************'''
|
|
1053
|
if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
|
|
|
1121
|
# x-axis point of the gaussian where the center is located from GaussFit of spectra
|
|
1054
|
try:
|
|
|
1122
|
GaussCenter = popt[1]
|
|
1055
|
popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
|
|
|
1123
|
ClosestCenter = xSamples_zoom[numpy.abs(xSamples_zoom-GaussCenter).argmin()]
|
|
1056
|
FitGauss=self.gaus(xSamples,*popt)
|
|
|
1124
|
PointGauCenter = numpy.where(xSamples_zoom==ClosestCenter)[0][0]
|
|
1057
|
|
|
|
|
|
|
1058
|
except :#RuntimeError:
|
|
|
|
|
|
1059
|
FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
|
|
|
|
|
|
1060
|
|
|
|
|
|
|
1061
|
|
|
|
|
|
|
1062
|
else:
|
|
|
|
|
|
1063
|
FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
|
|
|
1125
|
|
|
1064
|
|
|
|
1126
|
# Point where e^-1 is located in the gaussian
|
|
1065
|
|
|
|
|
|
|
1066
|
|
|
|
|
|
|
1067
|
'''****** Getting Fij ******'''
|
|
|
|
|
|
1068
|
Fijcspc = CSPCopt[:,2]/2*3
|
|
|
|
|
|
1069
|
|
|
|
|
|
|
1070
|
|
|
|
|
|
|
1071
|
GaussCenter = popt[1] #xFrec[GCpos]
|
|
|
|
|
|
1072
|
#Punto en Eje X de la Gaussiana donde se encuentra el centro
|
|
|
|
|
|
1073
|
ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
|
|
|
|
|
|
1074
|
PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
|
|
|
|
|
|
1075
|
|
|
|
|
|
|
1076
|
#Punto e^-1 hubicado en la Gaussiana
|
|
|
1127
|
PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
|
|
1077
|
PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
|
|
|
1128
|
FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # The closest point to"Peminus1" in "FitGauss"
|
|
1078
|
FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
|
|
|
1129
|
PointFij = numpy.where(FitGauss==FijClosest)[0][0]
|
|
1079
|
PointFij = numpy.where(FitGauss==FijClosest)[0][0]
|
|
|
1130
|
Fij = numpy.abs(xSamples_zoom[PointFij] - xSamples_zoom[PointGauCenter])
|
|
|
|
|
|
1131
|
|
|
1080
|
|
|
|
1132
|
'''********** Taking frequency ranges from mean SPCs **********'''
|
|
1081
|
if xSamples[PointFij] > xSamples[PointGauCenter]:
|
|
|
1133
|
GauWidth = popt[2] * 3/2 # Bandwidth of Gau01
|
|
1082
|
Fij = xSamples[PointFij] - xSamples[PointGauCenter]
|
|
|
|
|
|
1083
|
|
|
|
|
|
|
1084
|
else:
|
|
|
|
|
|
1085
|
Fij = xSamples[PointGauCenter] - xSamples[PointFij]
|
|
|
|
|
|
1086
|
|
|
|
|
|
|
1087
|
|
|
|
|
|
|
1088
|
'''****** Taking frequency ranges from SPCs ******'''
|
|
|
|
|
|
1089
|
|
|
|
|
|
|
1090
|
|
|
|
|
|
|
1091
|
#GaussCenter = popt[1] #Primer momento 01
|
|
|
|
|
|
1092
|
GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
|
|
|
1134
|
Range = numpy.empty(2)
|
|
1093
|
Range = numpy.empty(2)
|
|
|
1135
|
Range[0] = GaussCenter - GauWidth
|
|
1094
|
Range[0] = GaussCenter - GauWidth
|
|
|
1136
|
Range[1] = GaussCenter + GauWidth
|
|
1095
|
Range[1] = GaussCenter + GauWidth
|
|
|
1137
|
# Point in x-axis where the bandwidth is located (min:max)
|
|
1096
|
#Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
|
|
|
1138
|
ClosRangeMin = xSamples_zoom[numpy.abs(xSamples_zoom-Range[0]).argmin()]
|
|
1097
|
ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
|
|
|
1139
|
ClosRangeMax = xSamples_zoom[numpy.abs(xSamples_zoom-Range[1]).argmin()]
|
|
1098
|
ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
|
|
|
1140
|
PointRangeMin = numpy.where(xSamples_zoom==ClosRangeMin)[0][0]
|
|
1099
|
|
|
|
1141
|
PointRangeMax = numpy.where(xSamples_zoom==ClosRangeMax)[0][0]
|
|
1100
|
PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
|
|
|
|
|
|
1101
|
PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
|
|
|
|
|
|
1102
|
|
|
|
1142
|
Range=numpy.array([ PointRangeMin, PointRangeMax ])
|
|
1103
|
Range=numpy.array([ PointRangeMin, PointRangeMax ])
|
|
|
1143
|
FrecRange = xSamples_zoom[ Range[0] : Range[1] ]
|
|
|
|
|
|
1144
|
|
|
1104
|
|
|
|
1145
|
'''************************** Getting Phase Slope ***************************'''
|
|
1105
|
FrecRange = xFrec[ Range[0] : Range[1] ]
|
|
|
1146
|
for i in range(nPair):
|
|
1106
|
VelRange = xVel[ Range[0] : Range[1] ]
|
|
|
1147
|
if len(FrecRange) > 5:
|
|
1107
|
|
|
|
1148
|
PhaseRange = phase[i, xvalid[0][Range[0]:Range[1]]].copy()
|
|
1108
|
|
|
|
|
|
|
1109
|
'''****** Getting SCPC Slope ******'''
|
|
|
|
|
|
1110
|
|
|
|
|
|
|
1111
|
for i in range(spc.shape[0]):
|
|
|
|
|
|
1112
|
|
|
|
|
|
|
1113
|
if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
|
|
|
|
|
|
1114
|
PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
|
|
|
|
|
|
1115
|
|
|
|
|
|
|
1116
|
'''***********************VelRange******************'''
|
|
|
|
|
|
1117
|
|
|
|
1149
|
mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
|
|
1118
|
mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
|
|
|
|
|
|
1119
|
|
|
|
1150
|
if len(FrecRange) == len(PhaseRange):
|
|
1120
|
if len(FrecRange) == len(PhaseRange):
|
|
|
1151
|
try:
|
|
1121
|
try:
|
|
|
1152
|
slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
|
|
1122
|
slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
|
|
|
1153
|
PhaseSlope[i]=slope
|
|
1123
|
PhaseSlope[i]=slope
|
|
|
1154
|
PhaseInter[i]=intercept
|
|
1124
|
PhaseInter[i]=intercept
|
|
|
1155
|
except:
|
|
1125
|
except:
|
|
|
1156
|
return self.StopWindEstimation(error_code = 6)
|
|
1126
|
PhaseSlope[i]=0
|
|
|
|
|
|
1127
|
PhaseInter[i]=0
|
|
|
1157
|
else:
|
|
1128
|
else:
|
|
|
1158
|
return self.StopWindEstimation(error_code = 7)
|
|
1129
|
PhaseSlope[i]=0
|
|
|
|
|
|
1130
|
PhaseInter[i]=0
|
|
|
1159
|
else:
|
|
1131
|
else:
|
|
|
1160
|
return self.StopWindEstimation(error_code = 8)
|
|
1132
|
PhaseSlope[i]=0
|
|
|
|
|
|
1133
|
PhaseInter[i]=0
|
|
|
1161
|
|
|
1134
|
|
|
|
1162
|
'''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
|
|
|
|
|
|
1163
|
|
|
1135
|
|
|
|
1164
|
'''Getting constant C'''
|
|
1136
|
'''Getting constant C'''
|
|
|
1165
|
cC=(Fij*numpy.pi)**2
|
|
1137
|
cC=(Fij*numpy.pi)**2
|
|
|
1166
|
|
|
1138
|
|
|
|
1167
|
'''****** Getting constants F and G ******'''
|
|
1139
|
'''****** Getting constants F and G ******'''
|
|
|
1168
|
MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
|
|
1140
|
MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
|
|
|
1169
|
# MijEijNij = numpy.array([[Xi01,Eta01], [Xi02,Eta02], [Xi12,Eta12]])
|
|
1141
|
MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
|
|
|
1170
|
# MijResult0 = (-PhaseSlope[0] * cC) / (2*numpy.pi)
|
|
1142
|
MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
|
|
|
1171
|
MijResult1 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
|
|
1143
|
MijResults=numpy.array([MijResult0,MijResult1])
|
|
|
1172
|
MijResult2 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
|
|
|
|
|
|
1173
|
# MijResults = numpy.array([MijResult0, MijResult1, MijResult2])
|
|
|
|
|
|
1174
|
MijResults = numpy.array([MijResult1, MijResult2])
|
|
|
|
|
|
1175
|
(cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
|
|
1144
|
(cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
|
|
|
1176
|
|
|
1145
|
|
|
|
1177
|
'''****** Getting constants A, B and H ******'''
|
|
1146
|
'''****** Getting constants A, B and H ******'''
|
|
|
1178
|
W01 = numpy.nanmax( FitGauss01 )
|
|
1147
|
W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
|
|
|
1179
|
W02 = numpy.nanmax( FitGauss02 )
|
|
1148
|
W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
|
|
|
1180
|
W12 = numpy.nanmax( FitGauss12 )
|
|
1149
|
W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
|
|
|
|
|
|
1150
|
|
|
|
|
|
|
1151
|
WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
|
|
|
|
|
|
1152
|
WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
|
|
|
|
|
|
1153
|
WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
|
|
|
1181
|
|
|
1154
|
|
|
|
1182
|
WijResult01 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
|
|
1155
|
WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
|
|
|
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])
|
|
|
|
|
|
1186
|
|
|
1156
|
|
|
|
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] ])
|
|
1157
|
WijEijNij=numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
|
|
|
1188
|
(cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
|
|
1158
|
(cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
|
|
|
1189
|
|
|
1159
|
|
|
|
1190
|
VxVy=numpy.array([[cA,cH],[cH,cB]])
|
|
1160
|
VxVy=numpy.array([[cA,cH],[cH,cB]])
|
|
|
1191
|
VxVyResults=numpy.array([-cF,-cG])
|
|
1161
|
VxVyResults=numpy.array([-cF,-cG])
|
|
|
1192
|
(Vmer,Vzon) = numpy.linalg.solve(VxVy, VxVyResults)
|
|
1162
|
(Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
|
|
|
1193
|
Vver = -SPCMoments[1]*SPEED_OF_LIGHT/(2*radfreq)
|
|
1163
|
|
|
|
1194
|
error_code = 0
|
|
1164
|
Vzon = Vy
|
|
|
|
|
|
1165
|
Vmer = Vx
|
|
|
|
|
|
1166
|
Vmag=numpy.sqrt(Vzon**2+Vmer**2)
|
|
|
|
|
|
1167
|
Vang=numpy.arctan2(Vmer,Vzon)
|
|
|
|
|
|
1168
|
if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
|
|
|
|
|
|
1169
|
Vver=popt[1]
|
|
|
|
|
|
1170
|
else:
|
|
|
|
|
|
1171
|
Vver=numpy.NaN
|
|
|
|
|
|
1172
|
FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
|
|
|
|
|
|
1173
|
|
|
|
1195
|
|
|
1174
|
|
|
|
1196
|
return Vzon, Vmer, Vver, error_code
|
|
1175
|
return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
|
|
|
1197
|
|
|
1176
|
|
|
|
1198
|
class SpectralMoments(Operation):
|
|
1177
|
class SpectralMoments(Operation):
|
|
|
1199
|
|
|
1178
|
|
|
|
@@
-1220,12
+1199,13
class SpectralMoments(Operation):
|
|
|
1220
|
|
|
1199
|
|
|
|
1221
|
Affected:
|
|
1200
|
Affected:
|
|
|
1222
|
self.dataOut.moments : Parameters per channel
|
|
1201
|
self.dataOut.moments : Parameters per channel
|
|
|
1223
|
self.dataOut.data_snr : SNR per channel
|
|
1202
|
self.dataOut.data_SNR : SNR per channel
|
|
|
1224
|
|
|
1203
|
|
|
|
1225
|
'''
|
|
1204
|
'''
|
|
|
1226
|
|
|
1205
|
|
|
|
1227
|
def run(self, dataOut):
|
|
1206
|
def run(self, dataOut):
|
|
|
1228
|
|
|
1207
|
|
|
|
|
|
|
1208
|
#dataOut.data_pre = dataOut.data_pre[0]
|
|
|
1229
|
data = dataOut.data_pre[0]
|
|
1209
|
data = dataOut.data_pre[0]
|
|
|
1230
|
absc = dataOut.abscissaList[:-1]
|
|
1210
|
absc = dataOut.abscissaList[:-1]
|
|
|
1231
|
noise = dataOut.noise
|
|
1211
|
noise = dataOut.noise
|
|
|
@@
-1236,11
+1216,10
class SpectralMoments(Operation):
|
|
|
1236
|
data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
|
|
1216
|
data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
|
|
|
1237
|
|
|
1217
|
|
|
|
1238
|
dataOut.moments = data_param[:,1:,:]
|
|
1218
|
dataOut.moments = data_param[:,1:,:]
|
|
|
1239
|
dataOut.data_snr = data_param[:,0]
|
|
1219
|
dataOut.data_SNR = data_param[:,0]
|
|
|
1240
|
dataOut.data_pow = data_param[:,1]
|
|
1220
|
dataOut.data_POW = data_param[:,1]
|
|
|
1241
|
dataOut.data_dop = data_param[:,2]
|
|
1221
|
dataOut.data_DOP = data_param[:,2]
|
|
|
1242
|
dataOut.data_width = data_param[:,3]
|
|
1222
|
dataOut.data_WIDTH = data_param[:,3]
|
|
|
1243
|
|
|
|
|
|
|
1244
|
return dataOut
|
|
1223
|
return dataOut
|
|
|
1245
|
|
|
1224
|
|
|
|
1246
|
def __calculateMoments(self, oldspec, oldfreq, n0,
|
|
1225
|
def __calculateMoments(self, oldspec, oldfreq, n0,
|
|
|
@@
-1267,27
+1246,25
class SpectralMoments(Operation):
|
|
|
1267
|
vec_w = numpy.zeros(oldspec.shape[1])
|
|
1246
|
vec_w = numpy.zeros(oldspec.shape[1])
|
|
|
1268
|
vec_snr = numpy.zeros(oldspec.shape[1])
|
|
1247
|
vec_snr = numpy.zeros(oldspec.shape[1])
|
|
|
1269
|
|
|
1248
|
|
|
|
1270
|
# oldspec = numpy.ma.masked_invalid(oldspec)
|
|
1249
|
oldspec = numpy.ma.masked_invalid(oldspec)
|
|
|
1271
|
|
|
1250
|
|
|
|
1272
|
for ind in range(oldspec.shape[1]):
|
|
1251
|
for ind in range(oldspec.shape[1]):
|
|
|
1273
|
|
|
1252
|
|
|
|
1274
|
spec = oldspec[:,ind]
|
|
1253
|
spec = oldspec[:,ind]
|
|
|
1275
|
aux = spec*fwindow
|
|
1254
|
aux = spec*fwindow
|
|
|
1276
|
max_spec = aux.max()
|
|
1255
|
max_spec = aux.max()
|
|
|
1277
|
m = aux.tolist().index(max_spec)
|
|
1256
|
m = list(aux).index(max_spec)
|
|
|
1278
|
|
|
1257
|
|
|
|
1279
|
#Smooth
|
|
1258
|
#Smooth
|
|
|
1280
|
if (smooth == 0):
|
|
1259
|
if (smooth == 0): spec2 = spec
|
|
|
1281
|
spec2 = spec
|
|
1260
|
else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
|
|
|
1282
|
else:
|
|
|
|
|
|
1283
|
spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
|
|
|
|
|
|
1284
|
|
|
1261
|
|
|
|
1285
|
# Moments Estimation
|
|
1262
|
# Calculo de Momentos
|
|
|
1286
|
bb = spec2[numpy.arange(m,spec2.size)]
|
|
1263
|
bb = spec2[list(range(m,spec2.size))]
|
|
|
1287
|
bb = (bb<n0).nonzero()
|
|
1264
|
bb = (bb<n0).nonzero()
|
|
|
1288
|
bb = bb[0]
|
|
1265
|
bb = bb[0]
|
|
|
1289
|
|
|
1266
|
|
|
|
1290
|
ss = spec2[numpy.arange(0,m + 1)]
|
|
1267
|
ss = spec2[list(range(0,m + 1))]
|
|
|
1291
|
ss = (ss<n0).nonzero()
|
|
1268
|
ss = (ss<n0).nonzero()
|
|
|
1292
|
ss = ss[0]
|
|
1269
|
ss = ss[0]
|
|
|
1293
|
|
|
1270
|
|
|
|
@@
-1298,32
+1275,27
class SpectralMoments(Operation):
|
|
|
1298
|
if (bb0 < 0):
|
|
1275
|
if (bb0 < 0):
|
|
|
1299
|
bb0 = 0
|
|
1276
|
bb0 = 0
|
|
|
1300
|
|
|
1277
|
|
|
|
1301
|
if (ss.size == 0):
|
|
1278
|
if (ss.size == 0): ss1 = 1
|
|
|
1302
|
ss1 = 1
|
|
1279
|
else: ss1 = max(ss) + 1
|
|
|
1303
|
else:
|
|
|
|
|
|
1304
|
ss1 = max(ss) + 1
|
|
|
|
|
|
1305
|
|
|
|
|
|
|
1306
|
if (ss1 > m):
|
|
|
|
|
|
1307
|
ss1 = m
|
|
|
|
|
|
1308
|
|
|
1280
|
|
|
|
1309
|
valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
|
|
1281
|
if (ss1 > m): ss1 = m
|
|
|
1310
|
|
|
1282
|
|
|
|
1311
|
signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. Scipión added with correct definition
|
|
1283
|
valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1
|
|
|
1312
|
total_power = (spec2[valid] * fwindow[valid]).mean() # D. Scipión added with correct definition
|
|
|
|
|
|
1313
|
power = ((spec2[valid] - n0)*fwindow[valid]).sum()
|
|
1284
|
power = ((spec2[valid] - n0)*fwindow[valid]).sum()
|
|
|
1314
|
fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
|
|
1285
|
fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
|
|
|
1315
|
w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
|
|
1286
|
w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
|
|
|
1316
|
snr = (spec2.mean()-n0)/n0
|
|
1287
|
snr = (spec2.mean()-n0)/n0
|
|
|
|
|
|
1288
|
|
|
|
1317
|
if (snr < 1.e-20) :
|
|
1289
|
if (snr < 1.e-20) :
|
|
|
1318
|
snr = 1.e-20
|
|
1290
|
snr = 1.e-20
|
|
|
1319
|
|
|
1291
|
|
|
|
1320
|
# vec_power[ind] = power #D. Scipión replaced with the line below
|
|
1292
|
vec_power[ind] = power
|
|
|
1321
|
vec_power[ind] = total_power
|
|
|
|
|
|
1322
|
vec_fd[ind] = fd
|
|
1293
|
vec_fd[ind] = fd
|
|
|
1323
|
vec_w[ind] = w
|
|
1294
|
vec_w[ind] = w
|
|
|
1324
|
vec_snr[ind] = snr
|
|
1295
|
vec_snr[ind] = snr
|
|
|
1325
|
|
|
1296
|
|
|
|
1326
|
return numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
|
|
1297
|
moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
|
|
|
|
|
|
1298
|
return moments
|
|
|
1327
|
|
|
1299
|
|
|
|
1328
|
#------------------ Get SA Parameters --------------------------
|
|
1300
|
#------------------ Get SA Parameters --------------------------
|
|
|
1329
|
|
|
1301
|
|
|
|
@@
-1366,7
+1338,7
class SALags(Operation):
|
|
|
1366
|
self.dataOut.abscissaList
|
|
1338
|
self.dataOut.abscissaList
|
|
|
1367
|
self.dataOut.noise
|
|
1339
|
self.dataOut.noise
|
|
|
1368
|
self.dataOut.normFactor
|
|
1340
|
self.dataOut.normFactor
|
|
|
1369
|
self.dataOut.data_snr
|
|
1341
|
self.dataOut.data_SNR
|
|
|
1370
|
self.dataOut.groupList
|
|
1342
|
self.dataOut.groupList
|
|
|
1371
|
self.dataOut.nChannels
|
|
1343
|
self.dataOut.nChannels
|
|
|
1372
|
|
|
1344
|
|
|
|
@@
-1385,7
+1357,7
class SALags(Operation):
|
|
|
1385
|
nHeights = dataOut.nHeights
|
|
1357
|
nHeights = dataOut.nHeights
|
|
|
1386
|
absc = dataOut.abscissaList
|
|
1358
|
absc = dataOut.abscissaList
|
|
|
1387
|
noise = dataOut.noise
|
|
1359
|
noise = dataOut.noise
|
|
|
1388
|
SNR = dataOut.data_snr
|
|
1360
|
SNR = dataOut.data_SNR
|
|
|
1389
|
nChannels = dataOut.nChannels
|
|
1361
|
nChannels = dataOut.nChannels
|
|
|
1390
|
# pairsList = dataOut.groupList
|
|
1362
|
# pairsList = dataOut.groupList
|
|
|
1391
|
# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
|
|
1363
|
# pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
|
|
|
@@
-1455,6
+1427,11
class SALags(Operation):
|
|
|
1455
|
|
|
1427
|
|
|
|
1456
|
return phase
|
|
1428
|
return phase
|
|
|
1457
|
|
|
1429
|
|
|
|
|
|
|
1430
|
def fit_func( x, a0, a1, a2): #, a3, a4, a5):
|
|
|
|
|
|
1431
|
z = (x - a1) / a2
|
|
|
|
|
|
1432
|
y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
|
|
|
|
|
|
1433
|
return y
|
|
|
|
|
|
1434
|
|
|
|
1458
|
class SpectralFitting(Operation):
|
|
1435
|
class SpectralFitting(Operation):
|
|
|
1459
|
'''
|
|
1436
|
'''
|
|
|
1460
|
Function GetMoments()
|
|
1437
|
Function GetMoments()
|
|
|
@@
-1463,50
+1440,816
class SpectralFitting(Operation):
|
|
|
1463
|
Output:
|
|
1440
|
Output:
|
|
|
1464
|
Variables modified:
|
|
1441
|
Variables modified:
|
|
|
1465
|
'''
|
|
1442
|
'''
|
|
|
|
|
|
1443
|
def __calculateMoments(self,oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
|
|
|
1466
|
|
|
1444
|
|
|
|
1467
|
def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
|
|
1445
|
if (nicoh is None): nicoh = 1
|
|
|
|
|
|
1446
|
if (graph is None): graph = 0
|
|
|
|
|
|
1447
|
if (smooth is None): smooth = 0
|
|
|
|
|
|
1448
|
elif (self.smooth < 3): smooth = 0
|
|
|
|
|
|
1449
|
|
|
|
|
|
|
1450
|
if (type1 is None): type1 = 0
|
|
|
|
|
|
1451
|
if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
|
|
|
|
|
|
1452
|
if (snrth is None): snrth = -3
|
|
|
|
|
|
1453
|
if (dc is None): dc = 0
|
|
|
|
|
|
1454
|
if (aliasing is None): aliasing = 0
|
|
|
|
|
|
1455
|
if (oldfd is None): oldfd = 0
|
|
|
|
|
|
1456
|
if (wwauto is None): wwauto = 0
|
|
|
|
|
|
1457
|
|
|
|
|
|
|
1458
|
if (n0 < 1.e-20): n0 = 1.e-20
|
|
|
|
|
|
1459
|
|
|
|
|
|
|
1460
|
freq = oldfreq
|
|
|
|
|
|
1461
|
vec_power = numpy.zeros(oldspec.shape[1])
|
|
|
|
|
|
1462
|
vec_fd = numpy.zeros(oldspec.shape[1])
|
|
|
|
|
|
1463
|
vec_w = numpy.zeros(oldspec.shape[1])
|
|
|
|
|
|
1464
|
vec_snr = numpy.zeros(oldspec.shape[1])
|
|
|
|
|
|
1465
|
|
|
|
|
|
|
1466
|
oldspec = numpy.ma.masked_invalid(oldspec)
|
|
|
|
|
|
1467
|
|
|
|
|
|
|
1468
|
for ind in range(oldspec.shape[1]):
|
|
|
|
|
|
1469
|
|
|
|
|
|
|
1470
|
spec = oldspec[:,ind]
|
|
|
|
|
|
1471
|
aux = spec*fwindow
|
|
|
|
|
|
1472
|
max_spec = aux.max()
|
|
|
|
|
|
1473
|
m = list(aux).index(max_spec)
|
|
|
|
|
|
1474
|
|
|
|
|
|
|
1475
|
#Smooth
|
|
|
|
|
|
1476
|
if (smooth == 0): spec2 = spec
|
|
|
|
|
|
1477
|
else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
|
|
|
|
|
|
1478
|
|
|
|
|
|
|
1479
|
# Calculo de Momentos
|
|
|
|
|
|
1480
|
bb = spec2[list(range(m,spec2.size))]
|
|
|
|
|
|
1481
|
bb = (bb<n0).nonzero()
|
|
|
|
|
|
1482
|
bb = bb[0]
|
|
|
|
|
|
1483
|
|
|
|
|
|
|
1484
|
ss = spec2[list(range(0,m + 1))]
|
|
|
|
|
|
1485
|
ss = (ss<n0).nonzero()
|
|
|
|
|
|
1486
|
ss = ss[0]
|
|
|
|
|
|
1487
|
|
|
|
|
|
|
1488
|
if (bb.size == 0):
|
|
|
|
|
|
1489
|
bb0 = spec.size - 1 - m
|
|
|
|
|
|
1490
|
else:
|
|
|
|
|
|
1491
|
bb0 = bb[0] - 1
|
|
|
|
|
|
1492
|
if (bb0 < 0):
|
|
|
|
|
|
1493
|
bb0 = 0
|
|
|
|
|
|
1494
|
|
|
|
|
|
|
1495
|
if (ss.size == 0): ss1 = 1
|
|
|
|
|
|
1496
|
else: ss1 = max(ss) + 1
|
|
|
|
|
|
1497
|
|
|
|
|
|
|
1498
|
if (ss1 > m): ss1 = m
|
|
|
1468
|
|
|
1499
|
|
|
|
|
|
|
1500
|
valid = numpy.asarray(list(range(int(m + bb0 - ss1 + 1)))) + ss1
|
|
|
|
|
|
1501
|
power = ((spec2[valid] - n0)*fwindow[valid]).sum()
|
|
|
|
|
|
1502
|
fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
|
|
|
|
|
|
1503
|
w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
|
|
|
|
|
|
1504
|
snr = (spec2.mean()-n0)/n0
|
|
|
1469
|
|
|
1505
|
|
|
|
|
|
|
1506
|
if (snr < 1.e-20) :
|
|
|
|
|
|
1507
|
snr = 1.e-20
|
|
|
|
|
|
1508
|
|
|
|
|
|
|
1509
|
vec_power[ind] = power
|
|
|
|
|
|
1510
|
vec_fd[ind] = fd
|
|
|
|
|
|
1511
|
vec_w[ind] = w
|
|
|
|
|
|
1512
|
vec_snr[ind] = snr
|
|
|
|
|
|
1513
|
|
|
|
|
|
|
1514
|
moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
|
|
|
|
|
|
1515
|
return moments
|
|
|
|
|
|
1516
|
|
|
|
|
|
|
1517
|
#def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs):
|
|
|
|
|
|
1518
|
def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
|
|
|
|
|
|
1519
|
|
|
|
|
|
|
1520
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
1521
|
nProf = dataOut.nProfiles
|
|
|
|
|
|
1522
|
heights = dataOut.heightList
|
|
|
|
|
|
1523
|
nHei = len(heights)
|
|
|
|
|
|
1524
|
channels = dataOut.channelList
|
|
|
|
|
|
1525
|
nChan = len(channels)
|
|
|
|
|
|
1526
|
crosspairs = dataOut.groupList
|
|
|
|
|
|
1527
|
nPairs = len(crosspairs)
|
|
|
|
|
|
1528
|
#Separar espectros incoherentes de coherentes snr > 20 dB'
|
|
|
|
|
|
1529
|
snr_th = 10**(snrth/10.0)
|
|
|
|
|
|
1530
|
my_incoh_spectra = numpy.zeros([nChan, nProf,nHei], dtype='float')
|
|
|
|
|
|
1531
|
my_incoh_cspectra = numpy.zeros([nPairs,nProf, nHei], dtype='complex')
|
|
|
|
|
|
1532
|
my_incoh_aver = numpy.zeros([nChan, nHei])
|
|
|
|
|
|
1533
|
my_coh_aver = numpy.zeros([nChan, nHei])
|
|
|
|
|
|
1534
|
|
|
|
|
|
|
1535
|
coh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
|
|
|
|
|
|
1536
|
coh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
|
|
|
|
|
|
1537
|
coh_aver = numpy.zeros([nChan, nHei])
|
|
|
|
|
|
1538
|
|
|
|
|
|
|
1539
|
incoh_spectra = numpy.zeros([nChan, nProf, nHei], dtype='float')
|
|
|
|
|
|
1540
|
incoh_cspectra = numpy.zeros([nPairs, nProf, nHei], dtype='complex')
|
|
|
|
|
|
1541
|
incoh_aver = numpy.zeros([nChan, nHei])
|
|
|
|
|
|
1542
|
power = numpy.sum(spectra, axis=1)
|
|
|
|
|
|
1543
|
|
|
|
|
|
|
1544
|
if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
|
|
|
|
|
|
1545
|
if hei_th == None : hei_th = numpy.array([60,300,650])
|
|
|
|
|
|
1546
|
for ic in range(2):
|
|
|
|
|
|
1547
|
pair = crosspairs[ic]
|
|
|
|
|
|
1548
|
#si el SNR es mayor que el SNR threshold los datos se toman coherentes
|
|
|
|
|
|
1549
|
s_n0 = power[pair[0],:]/noise[pair[0]]
|
|
|
|
|
|
1550
|
s_n1 = power[pair[1],:]/noise[pair[1]]
|
|
|
|
|
|
1551
|
|
|
|
|
|
|
1552
|
valid1 =(s_n0>=snr_th).nonzero()
|
|
|
|
|
|
1553
|
valid2 = (s_n1>=snr_th).nonzero()
|
|
|
|
|
|
1554
|
#valid = valid2 + valid1 #numpy.concatenate((valid1,valid2), axis=None)
|
|
|
|
|
|
1555
|
valid1 = numpy.array(valid1[0])
|
|
|
|
|
|
1556
|
valid2 = numpy.array(valid2[0])
|
|
|
|
|
|
1557
|
valid = valid1
|
|
|
|
|
|
1558
|
for iv in range(len(valid2)):
|
|
|
|
|
|
1559
|
#for ivv in range(len(valid1)) :
|
|
|
|
|
|
1560
|
indv = numpy.array((valid1 == valid2[iv]).nonzero())
|
|
|
|
|
|
1561
|
if len(indv[0]) == 0 :
|
|
|
|
|
|
1562
|
valid = numpy.concatenate((valid,valid2[iv]), axis=None)
|
|
|
|
|
|
1563
|
if len(valid)>0:
|
|
|
|
|
|
1564
|
my_coh_aver[pair[0],valid]=1
|
|
|
|
|
|
1565
|
my_coh_aver[pair[1],valid]=1
|
|
|
|
|
|
1566
|
# si la coherencia es mayor a la coherencia threshold los datos se toman
|
|
|
|
|
|
1567
|
#print my_coh_aver[0,:]
|
|
|
|
|
|
1568
|
coh = numpy.squeeze(numpy.nansum(cspectra[ic,:,:], axis=0)/numpy.sqrt(numpy.nansum(spectra[pair[0],:,:], axis=0)*numpy.nansum(spectra[pair[1],:,:], axis=0)))
|
|
|
|
|
|
1569
|
#print('coh',numpy.absolute(coh))
|
|
|
|
|
|
1570
|
for ih in range(len(hei_th)):
|
|
|
|
|
|
1571
|
hvalid = (heights>hei_th[ih]).nonzero()
|
|
|
|
|
|
1572
|
hvalid = hvalid[0]
|
|
|
|
|
|
1573
|
if len(hvalid)>0:
|
|
|
|
|
|
1574
|
valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
|
|
|
|
|
|
1575
|
valid = valid[0]
|
|
|
|
|
|
1576
|
#print('hvalid:',hvalid)
|
|
|
|
|
|
1577
|
#print('valid', valid)
|
|
|
|
|
|
1578
|
if len(valid)>0:
|
|
|
|
|
|
1579
|
my_coh_aver[pair[0],hvalid[valid]] =1
|
|
|
|
|
|
1580
|
my_coh_aver[pair[1],hvalid[valid]] =1
|
|
|
|
|
|
1581
|
|
|
|
|
|
|
1582
|
coh_echoes = (my_coh_aver[pair[0],:] == 1).nonzero()
|
|
|
|
|
|
1583
|
incoh_echoes = (my_coh_aver[pair[0],:] != 1).nonzero()
|
|
|
|
|
|
1584
|
incoh_echoes = incoh_echoes[0]
|
|
|
|
|
|
1585
|
if len(incoh_echoes) > 0:
|
|
|
|
|
|
1586
|
my_incoh_spectra[pair[0],:,incoh_echoes] = spectra[pair[0],:,incoh_echoes]
|
|
|
|
|
|
1587
|
my_incoh_spectra[pair[1],:,incoh_echoes] = spectra[pair[1],:,incoh_echoes]
|
|
|
|
|
|
1588
|
my_incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
|
|
|
|
|
|
1589
|
my_incoh_aver[pair[0],incoh_echoes] = 1
|
|
|
|
|
|
1590
|
my_incoh_aver[pair[1],incoh_echoes] = 1
|
|
|
|
|
|
1591
|
|
|
|
|
|
|
1592
|
|
|
|
|
|
|
1593
|
for ic in range(2):
|
|
|
|
|
|
1594
|
pair = crosspairs[ic]
|
|
|
|
|
|
1595
|
|
|
|
|
|
|
1596
|
valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
|
|
|
|
|
|
1597
|
valid2 = (my_coh_aver[pair[1],:]==1).nonzero()
|
|
|
|
|
|
1598
|
valid1 = numpy.array(valid1[0])
|
|
|
|
|
|
1599
|
valid2 = numpy.array(valid2[0])
|
|
|
|
|
|
1600
|
valid = valid1
|
|
|
|
|
|
1601
|
#print valid1 , valid2
|
|
|
|
|
|
1602
|
for iv in range(len(valid2)):
|
|
|
|
|
|
1603
|
#for ivv in range(len(valid1)) :
|
|
|
|
|
|
1604
|
indv = numpy.array((valid1 == valid2[iv]).nonzero())
|
|
|
|
|
|
1605
|
if len(indv[0]) == 0 :
|
|
|
|
|
|
1606
|
valid = numpy.concatenate((valid,valid2[iv]), axis=None)
|
|
|
|
|
|
1607
|
#print valid
|
|
|
|
|
|
1608
|
#valid = numpy.concatenate((valid1,valid2), axis=None)
|
|
|
|
|
|
1609
|
valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
|
|
|
|
|
|
1610
|
valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
|
|
|
|
|
|
1611
|
valid1 = numpy.array(valid1[0])
|
|
|
|
|
|
1612
|
valid2 = numpy.array(valid2[0])
|
|
|
|
|
|
1613
|
incoh_echoes = valid1
|
|
|
|
|
|
1614
|
#print valid1, valid2
|
|
|
|
|
|
1615
|
#incoh_echoes= numpy.concatenate((valid1,valid2), axis=None)
|
|
|
|
|
|
1616
|
for iv in range(len(valid2)):
|
|
|
|
|
|
1617
|
#for ivv in range(len(valid1)) :
|
|
|
|
|
|
1618
|
indv = numpy.array((valid1 == valid2[iv]).nonzero())
|
|
|
|
|
|
1619
|
if len(indv[0]) == 0 :
|
|
|
|
|
|
1620
|
incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
|
|
|
|
|
|
1621
|
#print incoh_echoes
|
|
|
|
|
|
1622
|
if len(valid)>0:
|
|
|
|
|
|
1623
|
#print pair
|
|
|
|
|
|
1624
|
coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
|
|
|
|
|
|
1625
|
coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
|
|
|
|
|
|
1626
|
coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
|
|
|
|
|
|
1627
|
coh_aver[pair[0],valid]=1
|
|
|
|
|
|
1628
|
coh_aver[pair[1],valid]=1
|
|
|
|
|
|
1629
|
if len(incoh_echoes)>0:
|
|
|
|
|
|
1630
|
incoh_spectra[pair[0],:,incoh_echoes] = spectra[pair[0],:,incoh_echoes]
|
|
|
|
|
|
1631
|
incoh_spectra[pair[1],:,incoh_echoes] = spectra[pair[1],:,incoh_echoes]
|
|
|
|
|
|
1632
|
incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
|
|
|
|
|
|
1633
|
incoh_aver[pair[0],incoh_echoes]=1
|
|
|
|
|
|
1634
|
incoh_aver[pair[1],incoh_echoes]=1
|
|
|
|
|
|
1635
|
#plt.imshow(spectra[0,:,:],vmin=20000000)
|
|
|
|
|
|
1636
|
#plt.show()
|
|
|
|
|
|
1637
|
#my_incoh_aver = my_incoh_aver+1
|
|
|
|
|
|
1638
|
|
|
|
|
|
|
1639
|
#spec = my_incoh_spectra.copy()
|
|
|
|
|
|
1640
|
#cspec = my_incoh_cspectra.copy()
|
|
|
|
|
|
1641
|
#print('######################', spec)
|
|
|
|
|
|
1642
|
#print(self.numpy)
|
|
|
|
|
|
1643
|
#return spec, cspec,coh_aver
|
|
|
|
|
|
1644
|
return my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver
|
|
|
|
|
|
1645
|
|
|
|
|
|
|
1646
|
def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
|
|
|
|
|
|
1647
|
|
|
|
|
|
|
1648
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
1649
|
nProf = dataOut.nProfiles
|
|
|
|
|
|
1650
|
heights = dataOut.heightList
|
|
|
|
|
|
1651
|
nHei = len(heights)
|
|
|
|
|
|
1652
|
channels = dataOut.channelList
|
|
|
|
|
|
1653
|
nChan = len(channels)
|
|
|
|
|
|
1654
|
crosspairs = dataOut.groupList
|
|
|
|
|
|
1655
|
nPairs = len(crosspairs)
|
|
|
|
|
|
1656
|
|
|
|
|
|
|
1657
|
#data = dataOut.data_pre[0]
|
|
|
|
|
|
1658
|
absc = dataOut.abscissaList[:-1]
|
|
|
|
|
|
1659
|
#noise = dataOut.noise
|
|
|
|
|
|
1660
|
#nChannel = data.shape[0]
|
|
|
|
|
|
1661
|
data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
|
|
|
|
|
|
1662
|
|
|
|
|
|
|
1663
|
|
|
|
|
|
|
1664
|
#plt.plot(absc)
|
|
|
|
|
|
1665
|
#plt.show()
|
|
|
|
|
|
1666
|
clean_coh_spectra = spectra.copy()
|
|
|
|
|
|
1667
|
clean_coh_cspectra = cspectra.copy()
|
|
|
|
|
|
1668
|
clean_coh_aver = coh_aver.copy()
|
|
|
|
|
|
1669
|
|
|
|
|
|
|
1670
|
spwd_th=[10,6] #spwd_th[0] --> For satellites ; spwd_th[1] --> For special events like SUN.
|
|
|
|
|
|
1671
|
coh_th = 0.75
|
|
|
|
|
|
1672
|
|
|
|
|
|
|
1673
|
rtime0 = [6,18] # periodo sin ESF
|
|
|
|
|
|
1674
|
rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
|
|
|
|
|
|
1675
|
|
|
|
|
|
|
1676
|
time = index*5./60
|
|
|
|
|
|
1677
|
if clean_coh_echoes == 1 :
|
|
|
|
|
|
1678
|
for ind in range(nChan):
|
|
|
|
|
|
1679
|
data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
|
|
|
|
|
|
1680
|
#print data_param[:,3]
|
|
|
|
|
|
1681
|
spwd = data_param[:,3]
|
|
|
|
|
|
1682
|
#print spwd.shape
|
|
|
|
|
|
1683
|
# SPECB_JULIA,header=anal_header,jspectra=spectra,vel=velocities,hei=heights, num_aver=1, mode_fit=0,smoothing=smoothing,jvelr=velr,jspwd=spwd,jsnr=snr,jnoise=noise,jstdvnoise=stdvnoise
|
|
|
|
|
|
1684
|
#spwd1=[ 1.65607, 1.43416, 0.500373, 0.208361, 0.000000, 26.7767, 22.5936, 26.7530, 20.6962, 29.1098, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 28.0300, 27.0511, 27.8810, 26.3126, 27.8445, 24.6181, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000]
|
|
|
|
|
|
1685
|
#spwd=numpy.array([spwd1,spwd1,spwd1,spwd1])
|
|
|
|
|
|
1686
|
#print spwd.shape, heights.shape,coh_aver.shape
|
|
|
|
|
|
1687
|
# para obtener spwd
|
|
|
|
|
|
1688
|
for ic in range(nPairs):
|
|
|
|
|
|
1689
|
pair = crosspairs[ic]
|
|
|
|
|
|
1690
|
coh = numpy.squeeze(numpy.sum(cspectra[ic,:,:], axis=1)/numpy.sqrt(numpy.sum(spectra[pair[0],:,:], axis=1)*numpy.sum(spectra[pair[1],:,:], axis=1)))
|
|
|
|
|
|
1691
|
for ih in range(nHei) :
|
|
|
|
|
|
1692
|
# Considering heights higher than 200km in order to avoid removing phenomena like EEJ.
|
|
|
|
|
|
1693
|
if heights[ih] >= 200 and coh_aver[pair[0],ih] == 1 and coh_aver[pair[1],ih] == 1 :
|
|
|
|
|
|
1694
|
# Checking coherence
|
|
|
|
|
|
1695
|
if (numpy.abs(coh[ih]) <= coh_th) or (time >= rtime0[0] and time <= rtime0[1]) :
|
|
|
|
|
|
1696
|
# Checking spectral widths
|
|
|
|
|
|
1697
|
if (spwd[pair[0],ih] > spwd_th[0]) or (spwd[pair[1],ih] > spwd_th[0]) :
|
|
|
|
|
|
1698
|
# satelite
|
|
|
|
|
|
1699
|
clean_coh_spectra[pair,ih,:] = 0.0
|
|
|
|
|
|
1700
|
clean_coh_cspectra[ic,ih,:] = 0.0
|
|
|
|
|
|
1701
|
clean_coh_aver[pair,ih] = 0
|
|
|
|
|
|
1702
|
else :
|
|
|
|
|
|
1703
|
if ((spwd[pair[0],ih] < spwd_th[1]) or (spwd[pair[1],ih] < spwd_th[1])) :
|
|
|
|
|
|
1704
|
# Especial event like sun.
|
|
|
|
|
|
1705
|
clean_coh_spectra[pair,ih,:] = 0.0
|
|
|
|
|
|
1706
|
clean_coh_cspectra[ic,ih,:] = 0.0
|
|
|
|
|
|
1707
|
clean_coh_aver[pair,ih] = 0
|
|
|
|
|
|
1708
|
|
|
|
|
|
|
1709
|
return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
|
|
|
|
|
|
1710
|
|
|
|
|
|
|
1711
|
isConfig = False
|
|
|
|
|
|
1712
|
__dataReady = False
|
|
|
|
|
|
1713
|
bloques = None
|
|
|
|
|
|
1714
|
bloque0 = None
|
|
|
|
|
|
1715
|
|
|
|
|
|
|
1716
|
def __init__(self):
|
|
|
|
|
|
1717
|
Operation.__init__(self)
|
|
|
|
|
|
1718
|
self.i=0
|
|
|
|
|
|
1719
|
self.isConfig = False
|
|
|
|
|
|
1720
|
|
|
|
|
|
|
1721
|
|
|
|
|
|
|
1722
|
def setup(self,nChan,nProf,nHei,nBlocks):
|
|
|
|
|
|
1723
|
self.__dataReady = False
|
|
|
|
|
|
1724
|
self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
|
|
|
|
|
|
1725
|
self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
|
|
|
|
|
|
1726
|
|
|
|
|
|
|
1727
|
#def CleanRayleigh(self,dataOut,spectra,cspectra,out_spectra,out_cspectra,sat_spectra,sat_cspectra,crosspairs,heights, channels, nProf,nHei,nChan,nPairs,nIncohInt,nBlocks):
|
|
|
|
|
|
1728
|
def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
|
|
|
|
|
|
1729
|
#import matplotlib.pyplot as plt
|
|
|
|
|
|
1730
|
#for k in range(149):
|
|
|
|
|
|
1731
|
|
|
|
|
|
|
1732
|
# self.bloque0[:,:,:,k] = spectra[:,:,0:nHei]
|
|
|
|
|
|
1733
|
# self.bloques[:,:,:,k] = cspectra[:,:,0:nHei]
|
|
|
|
|
|
1734
|
#if self.i==nBlocks:
|
|
|
|
|
|
1735
|
# self.i==0
|
|
|
|
|
|
1736
|
rfunc = cspectra.copy() #self.bloques
|
|
|
|
|
|
1737
|
n_funct = len(rfunc[0,:,0,0])
|
|
|
|
|
|
1738
|
val_spc = spectra*0.0 #self.bloque0*0.0
|
|
|
|
|
|
1739
|
val_cspc = cspectra*0.0 #self.bloques*0.0
|
|
|
|
|
|
1740
|
in_sat_spectra = spectra.copy() #self.bloque0
|
|
|
|
|
|
1741
|
in_sat_cspectra = cspectra.copy() #self.bloques
|
|
|
|
|
|
1742
|
|
|
|
|
|
|
1743
|
#print( rfunc.shape)
|
|
|
|
|
|
1744
|
min_hei = 200
|
|
|
|
|
|
1745
|
nProf = dataOut.nProfiles
|
|
|
|
|
|
1746
|
heights = dataOut.heightList
|
|
|
|
|
|
1747
|
nHei = len(heights)
|
|
|
|
|
|
1748
|
channels = dataOut.channelList
|
|
|
|
|
|
1749
|
nChan = len(channels)
|
|
|
|
|
|
1750
|
crosspairs = dataOut.groupList
|
|
|
|
|
|
1751
|
nPairs = len(crosspairs)
|
|
|
|
|
|
1752
|
hval=(heights >= min_hei).nonzero()
|
|
|
|
|
|
1753
|
ih=hval[0]
|
|
|
|
|
|
1754
|
#print numpy.absolute(rfunc[:,0,0,14])
|
|
|
|
|
|
1755
|
for ih in range(hval[0][0],nHei):
|
|
|
|
|
|
1756
|
for ifreq in range(nProf):
|
|
|
|
|
|
1757
|
for ii in range(n_funct):
|
|
|
|
|
|
1758
|
|
|
|
|
|
|
1759
|
func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
|
|
|
|
|
|
1760
|
#print numpy.amin(func2clean)
|
|
|
|
|
|
1761
|
val = (numpy.isfinite(func2clean)==True).nonzero()
|
|
|
|
|
|
1762
|
if len(val)>0:
|
|
|
|
|
|
1763
|
min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
|
|
|
|
|
|
1764
|
if min_val <= -40 : min_val = -40
|
|
|
|
|
|
1765
|
max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
|
|
|
|
|
|
1766
|
if max_val >= 200 : max_val = 200
|
|
|
|
|
|
1767
|
#print min_val, max_val
|
|
|
|
|
|
1768
|
step = 1
|
|
|
|
|
|
1769
|
#Getting bins and the histogram
|
|
|
|
|
|
1770
|
x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
|
|
|
|
|
|
1771
|
y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
|
|
|
|
|
|
1772
|
mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
|
|
|
|
|
|
1773
|
sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
|
|
|
|
|
|
1774
|
parg = [numpy.amax(y_dist),mean,sigma]
|
|
|
|
|
|
1775
|
try :
|
|
|
|
|
|
1776
|
gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
|
|
|
|
|
|
1777
|
mode = gauss_fit[1]
|
|
|
|
|
|
1778
|
stdv = gauss_fit[2]
|
|
|
|
|
|
1779
|
except:
|
|
|
|
|
|
1780
|
mode = mean
|
|
|
|
|
|
1781
|
stdv = sigma
|
|
|
|
|
|
1782
|
# if ih == 14 and ii == 0 and ifreq ==0 :
|
|
|
|
|
|
1783
|
# print x_dist.shape, y_dist.shape
|
|
|
|
|
|
1784
|
# print x_dist, y_dist
|
|
|
|
|
|
1785
|
# print min_val, max_val, binstep
|
|
|
|
|
|
1786
|
# print func2clean
|
|
|
|
|
|
1787
|
# print mean,sigma
|
|
|
|
|
|
1788
|
# mean1,std = norm.fit(y_dist)
|
|
|
|
|
|
1789
|
# print mean1, std, gauss_fit
|
|
|
|
|
|
1790
|
# print fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
|
|
|
|
|
|
1791
|
# 7.84616 53.9307 3.61863
|
|
|
|
|
|
1792
|
#stdv = 3.61863 # 2.99089
|
|
|
|
|
|
1793
|
#mode = 53.9307 #7.79008
|
|
|
|
|
|
1794
|
|
|
|
|
|
|
1795
|
#Removing echoes greater than mode + 3*stdv
|
|
|
|
|
|
1796
|
factor_stdv = 2.5
|
|
|
|
|
|
1797
|
noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
|
|
|
|
|
|
1798
|
|
|
|
|
|
|
1799
|
if len(noval[0]) > 0:
|
|
|
|
|
|
1800
|
novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
|
|
|
|
|
|
1801
|
cross_pairs = crosspairs[ii]
|
|
|
|
|
|
1802
|
#Getting coherent echoes which are removed.
|
|
|
|
|
|
1803
|
if len(novall[0]) > 0:
|
|
|
|
|
|
1804
|
#val_spc[(0,1),novall[a],ih] = 1
|
|
|
|
|
|
1805
|
#val_spc[,(2,3),novall[a],ih] = 1
|
|
|
|
|
|
1806
|
val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
|
|
|
|
|
|
1807
|
val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
|
|
|
|
|
|
1808
|
val_cspc[novall[0],ii,ifreq,ih] = 1
|
|
|
|
|
|
1809
|
#print("OUT NOVALL 1")
|
|
|
|
|
|
1810
|
#Removing coherent from ISR data
|
|
|
|
|
|
1811
|
# if ih == 17 and ii == 0 and ifreq ==0 :
|
|
|
|
|
|
1812
|
# print spectra[:,cross_pairs[0],ifreq,ih]
|
|
|
|
|
|
1813
|
spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
|
|
|
|
|
|
1814
|
spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
|
|
|
|
|
|
1815
|
cspectra[noval,ii,ifreq,ih] = numpy.nan
|
|
|
|
|
|
1816
|
# if ih == 17 and ii == 0 and ifreq ==0 :
|
|
|
|
|
|
1817
|
# print spectra[:,cross_pairs[0],ifreq,ih]
|
|
|
|
|
|
1818
|
# print noval, len(noval[0])
|
|
|
|
|
|
1819
|
# print novall, len(novall[0])
|
|
|
|
|
|
1820
|
# print factor_stdv*stdv
|
|
|
|
|
|
1821
|
# print func2clean-mode
|
|
|
|
|
|
1822
|
# print val_spc[:,cross_pairs[0],ifreq,ih]
|
|
|
|
|
|
1823
|
# print spectra[:,cross_pairs[0],ifreq,ih]
|
|
|
|
|
|
1824
|
#no sale es para savedrifts >2
|
|
|
|
|
|
1825
|
''' channels = channels
|
|
|
|
|
|
1826
|
cross_pairs = cross_pairs
|
|
|
|
|
|
1827
|
#print("OUT NOVALL 2")
|
|
|
|
|
|
1828
|
|
|
|
|
|
|
1829
|
vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
|
|
|
|
|
|
1830
|
vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
|
|
|
|
|
|
1831
|
vcross = numpy.concatenate((vcross0,vcross1),axis=None)
|
|
|
|
|
|
1832
|
#print('vcros =', vcross)
|
|
|
|
|
|
1833
|
|
|
|
|
|
|
1834
|
#Getting coherent echoes which are removed.
|
|
|
|
|
|
1835
|
if len(novall) > 0:
|
|
|
|
|
|
1836
|
#val_spc[novall,ii,ifreq,ih] = 1
|
|
|
|
|
|
1837
|
val_spc[ii,ifreq,ih,novall] = 1
|
|
|
|
|
|
1838
|
if len(vcross) > 0:
|
|
|
|
|
|
1839
|
val_cspc[vcross,ifreq,ih,novall] = 1
|
|
|
|
|
|
1840
|
|
|
|
|
|
|
1841
|
#Removing coherent from ISR data.
|
|
|
|
|
|
1842
|
self.bloque0[ii,ifreq,ih,noval] = numpy.nan
|
|
|
|
|
|
1843
|
if len(vcross) > 0:
|
|
|
|
|
|
1844
|
self.bloques[vcross,ifreq,ih,noval] = numpy.nan
|
|
|
|
|
|
1845
|
'''
|
|
|
|
|
|
1846
|
#Getting average of the spectra and cross-spectra from incoherent echoes.
|
|
|
|
|
|
1847
|
out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
|
|
|
|
|
|
1848
|
out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
|
|
|
|
|
|
1849
|
for ih in range(nHei):
|
|
|
|
|
|
1850
|
for ifreq in range(nProf):
|
|
|
|
|
|
1851
|
for ich in range(nChan):
|
|
|
|
|
|
1852
|
tmp = spectra[:,ich,ifreq,ih]
|
|
|
|
|
|
1853
|
valid = (numpy.isfinite(tmp[:])==True).nonzero()
|
|
|
|
|
|
1854
|
# if ich == 0 and ifreq == 0 and ih == 17 :
|
|
|
|
|
|
1855
|
# print tmp
|
|
|
|
|
|
1856
|
# print valid
|
|
|
|
|
|
1857
|
# print len(valid[0])
|
|
|
|
|
|
1858
|
#print('TMP',tmp)
|
|
|
|
|
|
1859
|
if len(valid[0]) >0 :
|
|
|
|
|
|
1860
|
out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
|
|
|
|
|
|
1861
|
#for icr in range(nPairs):
|
|
|
|
|
|
1862
|
for icr in range(nPairs):
|
|
|
|
|
|
1863
|
tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
|
|
|
|
|
|
1864
|
valid = (numpy.isfinite(tmp)==True).nonzero()
|
|
|
|
|
|
1865
|
if len(valid[0]) > 0:
|
|
|
|
|
|
1866
|
out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
|
|
|
|
|
|
1867
|
# print('##########################################################')
|
|
|
|
|
|
1868
|
#Removing fake coherent echoes (at least 4 points around the point)
|
|
|
|
|
|
1869
|
|
|
|
|
|
|
1870
|
val_spectra = numpy.sum(val_spc,0)
|
|
|
|
|
|
1871
|
val_cspectra = numpy.sum(val_cspc,0)
|
|
|
|
|
|
1872
|
|
|
|
|
|
|
1873
|
val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
|
|
|
|
|
|
1874
|
val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
|
|
|
|
|
|
1875
|
|
|
|
|
|
|
1876
|
for i in range(nChan):
|
|
|
|
|
|
1877
|
for j in range(nProf):
|
|
|
|
|
|
1878
|
for k in range(nHei):
|
|
|
|
|
|
1879
|
if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
|
|
|
|
|
|
1880
|
val_spc[:,i,j,k] = 0.0
|
|
|
|
|
|
1881
|
for i in range(nPairs):
|
|
|
|
|
|
1882
|
for j in range(nProf):
|
|
|
|
|
|
1883
|
for k in range(nHei):
|
|
|
|
|
|
1884
|
if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
|
|
|
|
|
|
1885
|
val_cspc[:,i,j,k] = 0.0
|
|
|
|
|
|
1886
|
# val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHei*nChan))
|
|
|
|
|
|
1887
|
# if numpy.isfinite(val_spectra)==str(True):
|
|
|
|
|
|
1888
|
# noval = (val_spectra<1).nonzero()
|
|
|
|
|
|
1889
|
# if len(noval) > 0:
|
|
|
|
|
|
1890
|
# val_spc[:,noval] = 0.0
|
|
|
|
|
|
1891
|
# val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHei))
|
|
|
|
|
|
1892
|
|
|
|
|
|
|
1893
|
#val_cspc = numpy.reshape(val_spc, (149,nChan*nHei*nProf))
|
|
|
|
|
|
1894
|
#if numpy.isfinite(val_cspectra)==str(True):
|
|
|
|
|
|
1895
|
# noval = (val_cspectra<1).nonzero()
|
|
|
|
|
|
1896
|
# if len(noval) > 0:
|
|
|
|
|
|
1897
|
# val_cspc[:,noval] = 0.0
|
|
|
|
|
|
1898
|
# val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHei))
|
|
|
|
|
|
1899
|
|
|
|
|
|
|
1900
|
tmp_sat_spectra = spectra.copy()
|
|
|
|
|
|
1901
|
tmp_sat_spectra = tmp_sat_spectra*numpy.nan
|
|
|
|
|
|
1902
|
tmp_sat_cspectra = cspectra.copy()
|
|
|
|
|
|
1903
|
tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
|
|
|
|
|
|
1904
|
|
|
|
|
|
|
1905
|
# fig = plt.figure(figsize=(6,5))
|
|
|
|
|
|
1906
|
# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
|
|
|
|
|
|
1907
|
# ax = fig.add_axes([left, bottom, width, height])
|
|
|
|
|
|
1908
|
# cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
|
|
|
|
|
|
1909
|
# ax.clabel(cp, inline=True,fontsize=10)
|
|
|
|
|
|
1910
|
# plt.show()
|
|
|
|
|
|
1911
|
|
|
|
|
|
|
1912
|
val = (val_spc > 0).nonzero()
|
|
|
|
|
|
1913
|
if len(val[0]) > 0:
|
|
|
|
|
|
1914
|
tmp_sat_spectra[val] = in_sat_spectra[val]
|
|
|
|
|
|
1915
|
|
|
|
|
|
|
1916
|
val = (val_cspc > 0).nonzero()
|
|
|
|
|
|
1917
|
if len(val[0]) > 0:
|
|
|
|
|
|
1918
|
tmp_sat_cspectra[val] = in_sat_cspectra[val]
|
|
|
|
|
|
1919
|
|
|
|
|
|
|
1920
|
#Getting average of the spectra and cross-spectra from incoherent echoes.
|
|
|
|
|
|
1921
|
sat_spectra = numpy.zeros((nChan,nProf,nHei), dtype=float)
|
|
|
|
|
|
1922
|
sat_cspectra = numpy.zeros((nPairs,nProf,nHei), dtype=complex)
|
|
|
|
|
|
1923
|
for ih in range(nHei):
|
|
|
|
|
|
1924
|
for ifreq in range(nProf):
|
|
|
|
|
|
1925
|
for ich in range(nChan):
|
|
|
|
|
|
1926
|
tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
|
|
|
|
|
|
1927
|
valid = (numpy.isfinite(tmp)).nonzero()
|
|
|
|
|
|
1928
|
if len(valid[0]) > 0:
|
|
|
|
|
|
1929
|
sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
|
|
|
|
|
|
1930
|
|
|
|
|
|
|
1931
|
for icr in range(nPairs):
|
|
|
|
|
|
1932
|
tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
|
|
|
|
|
|
1933
|
valid = (numpy.isfinite(tmp)).nonzero()
|
|
|
|
|
|
1934
|
if len(valid[0]) > 0:
|
|
|
|
|
|
1935
|
sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
|
|
|
|
|
|
1936
|
#self.__dataReady= True
|
|
|
|
|
|
1937
|
#sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
|
|
|
|
|
|
1938
|
#if not self.__dataReady:
|
|
|
|
|
|
1939
|
#return None, None
|
|
|
|
|
|
1940
|
return out_spectra, out_cspectra,sat_spectra,sat_cspectra
|
|
|
|
|
|
1941
|
def REM_ISOLATED_POINTS(self,array,rth):
|
|
|
|
|
|
1942
|
# import matplotlib.pyplot as plt
|
|
|
|
|
|
1943
|
if rth == None : rth = 4
|
|
|
|
|
|
1944
|
|
|
|
|
|
|
1945
|
num_prof = len(array[0,:,0])
|
|
|
|
|
|
1946
|
num_hei = len(array[0,0,:])
|
|
|
|
|
|
1947
|
n2d = len(array[:,0,0])
|
|
|
|
|
|
1948
|
|
|
|
|
|
|
1949
|
for ii in range(n2d) :
|
|
|
|
|
|
1950
|
#print ii,n2d
|
|
|
|
|
|
1951
|
tmp = array[ii,:,:]
|
|
|
|
|
|
1952
|
#print tmp.shape, array[ii,101,:],array[ii,102,:]
|
|
|
|
|
|
1953
|
|
|
|
|
|
|
1954
|
# fig = plt.figure(figsize=(6,5))
|
|
|
|
|
|
1955
|
# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
|
|
|
|
|
|
1956
|
# ax = fig.add_axes([left, bottom, width, height])
|
|
|
|
|
|
1957
|
# x = range(num_prof)
|
|
|
|
|
|
1958
|
# y = range(num_hei)
|
|
|
|
|
|
1959
|
# cp = ax.contour(y,x,tmp)
|
|
|
|
|
|
1960
|
# ax.clabel(cp, inline=True,fontsize=10)
|
|
|
|
|
|
1961
|
# plt.show()
|
|
|
|
|
|
1962
|
|
|
|
|
|
|
1963
|
#indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
|
|
|
|
|
|
1964
|
tmp = numpy.reshape(tmp,num_prof*num_hei)
|
|
|
|
|
|
1965
|
indxs1 = (numpy.isfinite(tmp)==True).nonzero()
|
|
|
|
|
|
1966
|
indxs2 = (tmp > 0).nonzero()
|
|
|
|
|
|
1967
|
|
|
|
|
|
|
1968
|
indxs1 = (indxs1[0])
|
|
|
|
|
|
1969
|
indxs2 = indxs2[0]
|
|
|
|
|
|
1970
|
#indxs1 = numpy.array(indxs1[0])
|
|
|
|
|
|
1971
|
#indxs2 = numpy.array(indxs2[0])
|
|
|
|
|
|
1972
|
indxs = None
|
|
|
|
|
|
1973
|
#print indxs1 , indxs2
|
|
|
|
|
|
1974
|
for iv in range(len(indxs2)):
|
|
|
|
|
|
1975
|
indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
|
|
|
|
|
|
1976
|
#print len(indxs2), indv
|
|
|
|
|
|
1977
|
if len(indv[0]) > 0 :
|
|
|
|
|
|
1978
|
indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
|
|
|
|
|
|
1979
|
# print indxs
|
|
|
|
|
|
1980
|
indxs = indxs[1:]
|
|
|
|
|
|
1981
|
#print indxs, len(indxs)
|
|
|
|
|
|
1982
|
if len(indxs) < 4 :
|
|
|
|
|
|
1983
|
array[ii,:,:] = 0.
|
|
|
|
|
|
1984
|
return
|
|
|
|
|
|
1985
|
|
|
|
|
|
|
1986
|
xpos = numpy.mod(indxs ,num_hei)
|
|
|
|
|
|
1987
|
ypos = (indxs / num_hei)
|
|
|
|
|
|
1988
|
sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
|
|
|
|
|
|
1989
|
#print sx
|
|
|
|
|
|
1990
|
xpos = xpos[sx]
|
|
|
|
|
|
1991
|
ypos = ypos[sx]
|
|
|
|
|
|
1992
|
|
|
|
|
|
|
1993
|
# *********************************** Cleaning isolated points **********************************
|
|
|
|
|
|
1994
|
ic = 0
|
|
|
|
|
|
1995
|
while True :
|
|
|
|
|
|
1996
|
r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
|
|
|
|
|
|
1997
|
#no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
|
|
|
|
|
|
1998
|
#plt.plot(r)
|
|
|
|
|
|
1999
|
#plt.show()
|
|
|
|
|
|
2000
|
no_coh1 = (numpy.isfinite(r)==True).nonzero()
|
|
|
|
|
|
2001
|
no_coh2 = (r <= rth).nonzero()
|
|
|
|
|
|
2002
|
#print r, no_coh1, no_coh2
|
|
|
|
|
|
2003
|
no_coh1 = numpy.array(no_coh1[0])
|
|
|
|
|
|
2004
|
no_coh2 = numpy.array(no_coh2[0])
|
|
|
|
|
|
2005
|
no_coh = None
|
|
|
|
|
|
2006
|
#print valid1 , valid2
|
|
|
|
|
|
2007
|
for iv in range(len(no_coh2)):
|
|
|
|
|
|
2008
|
indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
|
|
|
|
|
|
2009
|
if len(indv[0]) > 0 :
|
|
|
|
|
|
2010
|
no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
|
|
|
|
|
|
2011
|
no_coh = no_coh[1:]
|
|
|
|
|
|
2012
|
#print len(no_coh), no_coh
|
|
|
|
|
|
2013
|
if len(no_coh) < 4 :
|
|
|
|
|
|
2014
|
#print xpos[ic], ypos[ic], ic
|
|
|
|
|
|
2015
|
# plt.plot(r)
|
|
|
|
|
|
2016
|
# plt.show()
|
|
|
|
|
|
2017
|
xpos[ic] = numpy.nan
|
|
|
|
|
|
2018
|
ypos[ic] = numpy.nan
|
|
|
|
|
|
2019
|
|
|
|
|
|
|
2020
|
ic = ic + 1
|
|
|
|
|
|
2021
|
if (ic == len(indxs)) :
|
|
|
|
|
|
2022
|
break
|
|
|
|
|
|
2023
|
#print( xpos, ypos)
|
|
|
|
|
|
2024
|
|
|
|
|
|
|
2025
|
indxs = (numpy.isfinite(list(xpos))==True).nonzero()
|
|
|
|
|
|
2026
|
#print indxs[0]
|
|
|
|
|
|
2027
|
if len(indxs[0]) < 4 :
|
|
|
|
|
|
2028
|
array[ii,:,:] = 0.
|
|
|
|
|
|
2029
|
return
|
|
|
|
|
|
2030
|
|
|
|
|
|
|
2031
|
xpos = xpos[indxs[0]]
|
|
|
|
|
|
2032
|
ypos = ypos[indxs[0]]
|
|
|
|
|
|
2033
|
for i in range(0,len(ypos)):
|
|
|
|
|
|
2034
|
ypos[i]=int(ypos[i])
|
|
|
|
|
|
2035
|
junk = tmp
|
|
|
|
|
|
2036
|
tmp = junk*0.0
|
|
|
|
|
|
2037
|
|
|
|
|
|
|
2038
|
tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
|
|
|
|
|
|
2039
|
array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
|
|
|
|
|
|
2040
|
|
|
|
|
|
|
2041
|
#print array.shape
|
|
|
|
|
|
2042
|
#tmp = numpy.reshape(tmp,(num_prof,num_hei))
|
|
|
|
|
|
2043
|
#print tmp.shape
|
|
|
|
|
|
2044
|
|
|
|
|
|
|
2045
|
# fig = plt.figure(figsize=(6,5))
|
|
|
|
|
|
2046
|
# left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
|
|
|
|
|
|
2047
|
# ax = fig.add_axes([left, bottom, width, height])
|
|
|
|
|
|
2048
|
# x = range(num_prof)
|
|
|
|
|
|
2049
|
# y = range(num_hei)
|
|
|
|
|
|
2050
|
# cp = ax.contour(y,x,array[ii,:,:])
|
|
|
|
|
|
2051
|
# ax.clabel(cp, inline=True,fontsize=10)
|
|
|
|
|
|
2052
|
# plt.show()
|
|
|
|
|
|
2053
|
return array
|
|
|
|
|
|
2054
|
def moments(self,doppler,yarray,npoints):
|
|
|
|
|
|
2055
|
ytemp = yarray
|
|
|
|
|
|
2056
|
#val = WHERE(ytemp GT 0,cval)
|
|
|
|
|
|
2057
|
#if cval == 0 : val = range(npoints-1)
|
|
|
|
|
|
2058
|
val = (ytemp > 0).nonzero()
|
|
|
|
|
|
2059
|
val = val[0]
|
|
|
|
|
|
2060
|
#print('hvalid:',hvalid)
|
|
|
|
|
|
2061
|
#print('valid', valid)
|
|
|
|
|
|
2062
|
if len(val) == 0 : val = range(npoints-1)
|
|
|
|
|
|
2063
|
|
|
|
|
|
|
2064
|
ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
|
|
|
|
|
|
2065
|
ytemp[len(ytemp):] = [ynew]
|
|
|
|
|
|
2066
|
|
|
|
|
|
|
2067
|
index = 0
|
|
|
|
|
|
2068
|
index = numpy.argmax(ytemp)
|
|
|
|
|
|
2069
|
ytemp = numpy.roll(ytemp,int(npoints/2)-1-index)
|
|
|
|
|
|
2070
|
ytemp = ytemp[0:npoints-1]
|
|
|
|
|
|
2071
|
|
|
|
|
|
|
2072
|
fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
|
|
|
|
|
|
2073
|
smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
|
|
|
|
|
|
2074
|
return [fmom,numpy.sqrt(smom)]
|
|
|
|
|
|
2075
|
# **********************************************************************************************
|
|
|
|
|
|
2076
|
index = 0
|
|
|
|
|
|
2077
|
fint = 0
|
|
|
|
|
|
2078
|
buffer = 0
|
|
|
|
|
|
2079
|
buffer2 = 0
|
|
|
|
|
|
2080
|
buffer3 = 0
|
|
|
|
|
|
2081
|
def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
|
|
|
|
|
|
2082
|
#print (dataOut.utctime)
|
|
|
|
|
|
2083
|
import matplotlib.pyplot as plt
|
|
|
|
|
|
2084
|
#nGroups = groupArray.shape[0]
|
|
|
|
|
|
2085
|
nChannels = dataOut.nChannels
|
|
|
|
|
|
2086
|
nHeights= dataOut.heightList.size
|
|
|
|
|
|
2087
|
nProf = dataOut.nProfiles
|
|
|
|
|
|
2088
|
|
|
|
|
|
|
2089
|
tini=time.localtime(dataOut.utctime)
|
|
|
|
|
|
2090
|
if (tini.tm_min % 5) == 0 and (tini.tm_sec < 5 and self.fint==0):
|
|
|
|
|
|
2091
|
# print tini.tm_min
|
|
|
|
|
|
2092
|
self.index = 0
|
|
|
|
|
|
2093
|
jspc = self.buffer
|
|
|
|
|
|
2094
|
jcspc = self.buffer2
|
|
|
|
|
|
2095
|
jnoise = self.buffer3
|
|
|
|
|
|
2096
|
|
|
|
|
|
|
2097
|
self.buffer = dataOut.data_spc
|
|
|
|
|
|
2098
|
self.buffer2 = dataOut.data_cspc
|
|
|
|
|
|
2099
|
self.buffer3 = dataOut.noise
|
|
|
|
|
|
2100
|
self.fint = 1
|
|
|
|
|
|
2101
|
#print self.buffer[0,:,0]
|
|
|
|
|
|
2102
|
|
|
|
|
|
|
2103
|
if numpy.any(jspc) :
|
|
|
|
|
|
2104
|
#print (len(jspc), jspc.shape)
|
|
|
|
|
|
2105
|
#print jspc[len(jspc)-4,:,0]
|
|
|
|
|
|
2106
|
jspc= numpy.reshape(jspc,(int(len(jspc)/4),nChannels,nProf,nHeights))
|
|
|
|
|
|
2107
|
jcspc= numpy.reshape(jcspc,(int(len(jcspc)/2),2,nProf,nHeights))
|
|
|
|
|
|
2108
|
jnoise= numpy.reshape(jnoise,(int(len(jnoise)/4),nChannels))
|
|
|
|
|
|
2109
|
#print jspc[len(jspc)-1,0,:,0]
|
|
|
|
|
|
2110
|
else:
|
|
|
|
|
|
2111
|
dataOut.flagNoData = True
|
|
|
|
|
|
2112
|
return dataOut
|
|
|
|
|
|
2113
|
|
|
|
|
|
|
2114
|
else :
|
|
|
|
|
|
2115
|
#print tini.tm_min
|
|
|
|
|
|
2116
|
#self.fint = 0
|
|
|
|
|
|
2117
|
if (tini.tm_min % 5) == 0 : self.fint = 1
|
|
|
|
|
|
2118
|
else : self.fint = 0
|
|
|
|
|
|
2119
|
self.index += 1
|
|
|
|
|
|
2120
|
#print( len(self.buffer))
|
|
|
|
|
|
2121
|
|
|
|
|
|
|
2122
|
if numpy.any(self.buffer):
|
|
|
|
|
|
2123
|
self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
|
|
|
|
|
|
2124
|
self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
|
|
|
|
|
|
2125
|
self.buffer3 = numpy.concatenate((self.buffer3,dataOut.noise), axis=0)
|
|
|
|
|
|
2126
|
else:
|
|
|
|
|
|
2127
|
self.buffer = dataOut.data_spc
|
|
|
|
|
|
2128
|
self.buffer2 = dataOut.data_cspc
|
|
|
|
|
|
2129
|
self.buffer3 = dataOut.noise
|
|
|
|
|
|
2130
|
#print self.index, self.fint
|
|
|
|
|
|
2131
|
#print self.buffer2.shape
|
|
|
|
|
|
2132
|
dataOut.flagNoData = True
|
|
|
|
|
|
2133
|
return dataOut
|
|
|
|
|
|
2134
|
# if self.index == 0 and self.fint == 1 :
|
|
|
|
|
|
2135
|
# if jspc != None:
|
|
|
|
|
|
2136
|
# print len(jspc), jspc.shape
|
|
|
|
|
|
2137
|
# jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4))
|
|
|
|
|
|
2138
|
# print jspc.shape
|
|
|
|
|
|
2139
|
# dataOut.flagNoData = True
|
|
|
|
|
|
2140
|
# return dataOut
|
|
|
1470
|
if path != None:
|
|
2141
|
if path != None:
|
|
|
1471
|
sys.path.append(path)
|
|
2142
|
sys.path.append(path)
|
|
|
1472
|
self.dataOut.library = importlib.import_module(file)
|
|
2143
|
self.library = importlib.import_module(file)
|
|
|
1473
|
|
|
2144
|
|
|
|
1474
|
#To be inserted as a parameter
|
|
2145
|
#To be inserted as a parameter
|
|
|
1475
|
groupArray = numpy.array(groupList)
|
|
2146
|
groupArray = numpy.array(groupList)
|
|
|
1476
|
#groupArray = numpy.array([[0,1],[2,3]])
|
|
2147
|
#groupArray = numpy.array([[0,1],[2,3]])
|
|
|
1477
|
self.dataOut.groupList = groupArray
|
|
2148
|
dataOut.groupList = groupArray
|
|
|
1478
|
|
|
2149
|
|
|
|
1479
|
nGroups = groupArray.shape[0]
|
|
2150
|
nGroups = groupArray.shape[0]
|
|
|
1480
|
nChannels = self.dataIn.nChannels
|
|
2151
|
nChannels = dataOut.nChannels
|
|
|
1481
|
nHeights=self.dataIn.heightList.size
|
|
2152
|
nHeights= dataOut.heightList.size
|
|
|
1482
|
|
|
2153
|
# print self.index
|
|
|
1483
|
#Parameters Array
|
|
2154
|
#Parameters Array
|
|
|
1484
|
self.dataOut.data_param = None
|
|
2155
|
dataOut.data_param = None
|
|
|
|
|
|
2156
|
dataOut.data_paramC = None
|
|
|
1485
|
|
|
2157
|
|
|
|
1486
|
#Set constants
|
|
2158
|
#Set constants
|
|
|
1487
|
constants = self.dataOut.library.setConstants(self.dataIn)
|
|
2159
|
constants = self.library.setConstants(dataOut)
|
|
|
1488
|
self.dataOut.constants = constants
|
|
2160
|
dataOut.constants = constants
|
|
|
1489
|
M = self.dataIn.normFactor
|
|
2161
|
M = dataOut.normFactor
|
|
|
1490
|
N = self.dataIn.nFFTPoints
|
|
2162
|
N = dataOut.nFFTPoints
|
|
|
1491
|
ippSeconds = self.dataIn.ippSeconds
|
|
2163
|
ippSeconds = dataOut.ippSeconds
|
|
|
1492
|
K = self.dataIn.nIncohInt
|
|
2164
|
K = dataOut.nIncohInt
|
|
|
1493
|
pairsArray = numpy.array(self.dataIn.pairsList)
|
|
2165
|
pairsArray = numpy.array(dataOut.pairsList)
|
|
|
1494
|
|
|
2166
|
|
|
|
|
|
|
2167
|
snrth= 20
|
|
|
|
|
|
2168
|
spectra = dataOut.data_spc
|
|
|
|
|
|
2169
|
cspectra = dataOut.data_cspc
|
|
|
|
|
|
2170
|
nProf = dataOut.nProfiles
|
|
|
|
|
|
2171
|
heights = dataOut.heightList
|
|
|
|
|
|
2172
|
nHei = len(heights)
|
|
|
|
|
|
2173
|
|
|
|
|
|
|
2174
|
channels = dataOut.channelList
|
|
|
|
|
|
2175
|
nChan = len(channels)
|
|
|
|
|
|
2176
|
nIncohInt = dataOut.nIncohInt
|
|
|
|
|
|
2177
|
crosspairs = dataOut.groupList
|
|
|
|
|
|
2178
|
noise = dataOut.noise
|
|
|
|
|
|
2179
|
#print( nProf,heights)
|
|
|
|
|
|
2180
|
#print( jspc.shape, jspc.shape[0])
|
|
|
|
|
|
2181
|
#print noise
|
|
|
|
|
|
2182
|
#print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise)
|
|
|
|
|
|
2183
|
jnoise = jnoise/N
|
|
|
|
|
|
2184
|
noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
|
|
|
|
|
|
2185
|
#print( noise)
|
|
|
|
|
|
2186
|
power = numpy.sum(spectra, axis=1)
|
|
|
|
|
|
2187
|
#print power[0,:]
|
|
|
|
|
|
2188
|
#print("CROSSPAIRS",crosspairs)
|
|
|
|
|
|
2189
|
nPairs = len(crosspairs)
|
|
|
|
|
|
2190
|
#print(numpy.shape(dataOut.data_spc))
|
|
|
|
|
|
2191
|
absc = dataOut.abscissaList[:-1]
|
|
|
|
|
|
2192
|
#print absc.shape
|
|
|
|
|
|
2193
|
#nBlocks=149
|
|
|
|
|
|
2194
|
#print('spectra', spectra.shape)
|
|
|
|
|
|
2195
|
#print('noise print', crosspairs)
|
|
|
|
|
|
2196
|
#print('spectra', spectra.shape)
|
|
|
|
|
|
2197
|
#print('cspectra', cspectra.shape)
|
|
|
|
|
|
2198
|
#print numpy.array(dataOut.data_pre[1]).shape
|
|
|
|
|
|
2199
|
#spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs)
|
|
|
|
|
|
2200
|
|
|
|
|
|
|
2201
|
if not self.isConfig:
|
|
|
|
|
|
2202
|
#self.setup(nChan,nProf,nHei=35,nBlocks=nBlocks)
|
|
|
|
|
|
2203
|
self.isConfig = True
|
|
|
|
|
|
2204
|
|
|
|
|
|
|
2205
|
#print ("configure todo")
|
|
|
|
|
|
2206
|
# dataOut.flagNoData = True
|
|
|
|
|
|
2207
|
index = tini.tm_hour*12+tini.tm_min/5
|
|
|
|
|
|
2208
|
#print index
|
|
|
|
|
|
2209
|
jspc = jspc/N/N
|
|
|
|
|
|
2210
|
jcspc = jcspc/N/N
|
|
|
|
|
|
2211
|
#dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
|
|
|
|
|
|
2212
|
tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
|
|
|
|
|
|
2213
|
jspectra = tmp_spectra*len(jspc[:,0,0,0])
|
|
|
|
|
|
2214
|
jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
|
|
|
|
|
|
2215
|
#incoh_spectra, incoh_cspectra,coh_aver = self.__DiffCoherent(snrth, dataOut.data_spc, dataOut.data_cspc, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs)
|
|
|
|
|
|
2216
|
my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth, None, None)
|
|
|
|
|
|
2217
|
clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
|
|
|
|
|
|
2218
|
dataOut.data_spc = incoh_spectra
|
|
|
|
|
|
2219
|
dataOut.data_cspc = incoh_cspectra
|
|
|
|
|
|
2220
|
#dataOut.data_spc = tmp_spectra
|
|
|
|
|
|
2221
|
#dataOut.data_cspc = tmp_cspectra
|
|
|
|
|
|
2222
|
|
|
|
|
|
|
2223
|
clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
|
|
|
|
|
|
2224
|
coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
|
|
|
|
|
|
2225
|
#plt.plot( tmp_spectra[0,:,17])
|
|
|
|
|
|
2226
|
#plt.show()
|
|
|
|
|
|
2227
|
# plt.plot( incoh_spectra[0,64,:])
|
|
|
|
|
|
2228
|
# plt.show()
|
|
|
|
|
|
2229
|
|
|
|
|
|
|
2230
|
# plt.imshow(dataOut.data_spc[0,:,:],vmin=20000000)
|
|
|
|
|
|
2231
|
# plt.show()
|
|
|
1495
|
#List of possible combinations
|
|
2232
|
#List of possible combinations
|
|
|
1496
|
listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
|
|
2233
|
listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
|
|
|
1497
|
indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
|
|
2234
|
indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
|
|
|
|
|
|
2235
|
#print("listComb",listComb)
|
|
|
1498
|
|
|
2236
|
|
|
|
1499
|
if getSNR:
|
|
2237
|
if getSNR:
|
|
|
1500
|
listChannels = groupArray.reshape((groupArray.size))
|
|
2238
|
listChannels = groupArray.reshape((groupArray.size))
|
|
|
1501
|
listChannels.sort()
|
|
2239
|
listChannels.sort()
|
|
|
1502
|
noise = self.dataIn.getNoise()
|
|
2240
|
#noise = dataOut.getNoise()
|
|
|
1503
|
self.dataOut.data_snr = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
|
|
2241
|
#print noise
|
|
|
1504
|
|
|
2242
|
#print(numpy.shape(noise))
|
|
|
|
|
|
2243
|
#dataOut.data_spc, dataOut.data_cspc = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights, nHei, nChan, channels, noise, nPairs, crosspairs)
|
|
|
|
|
|
2244
|
dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
|
|
|
|
|
|
2245
|
#dataOut.data_SNR = self.__getSNR(jspectra[listChannels,:,:], noise[listChannels])
|
|
|
|
|
|
2246
|
|
|
|
|
|
|
2247
|
if dataOut.data_paramC is None:
|
|
|
|
|
|
2248
|
dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
|
|
|
1505
|
for i in range(nGroups):
|
|
2249
|
for i in range(nGroups):
|
|
|
1506
|
coord = groupArray[i,:]
|
|
2250
|
coord = groupArray[i,:]
|
|
|
1507
|
|
|
|
|
|
|
1508
|
#Input data array
|
|
2251
|
#Input data array
|
|
|
1509
|
data = self.dataIn.data_spc[coord,:,:]/(M*N)
|
|
2252
|
data = dataOut.data_spc[coord,:,:]/(M*N)
|
|
|
1510
|
data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
|
|
2253
|
data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
|
|
|
1511
|
|
|
2254
|
|
|
|
1512
|
#Cross Spectra data array for Covariance Matrixes
|
|
2255
|
#Cross Spectra data array for Covariance Matrixes
|
|
|
@@
-1515,16
+2258,39
class SpectralFitting(Operation):
|
|
|
1515
|
pairsSel = numpy.array([coord[x],coord[y]])
|
|
2258
|
pairsSel = numpy.array([coord[x],coord[y]])
|
|
|
1516
|
indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
|
|
2259
|
indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
|
|
|
1517
|
ind += 1
|
|
2260
|
ind += 1
|
|
|
1518
|
dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
|
|
2261
|
dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
|
|
|
1519
|
dataCross = dataCross**2/K
|
|
2262
|
dataCross = dataCross**2
|
|
|
|
|
|
2263
|
#print dataOut.data_SNR.shape
|
|
|
|
|
|
2264
|
|
|
|
|
|
|
2265
|
nhei = nHeights
|
|
|
|
|
|
2266
|
poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
|
|
|
|
|
|
2267
|
if i == 0 : my_noises = numpy.zeros(4,dtype=float) #FLTARR(4)
|
|
|
|
|
|
2268
|
n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
|
|
|
|
|
|
2269
|
n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
|
|
|
|
|
|
2270
|
|
|
|
|
|
|
2271
|
n0 = n0i
|
|
|
|
|
|
2272
|
n1= n1i
|
|
|
|
|
|
2273
|
my_noises[2*i+0] = n0
|
|
|
|
|
|
2274
|
my_noises[2*i+1] = n1
|
|
|
|
|
|
2275
|
snrth = -16.0
|
|
|
|
|
|
2276
|
snrth = 10**(snrth/10.0)
|
|
|
1520
|
|
|
2277
|
|
|
|
1521
|
for h in range(nHeights):
|
|
2278
|
for h in range(nHeights):
|
|
|
1522
|
|
|
2279
|
# print("I ", "H", i,h )
|
|
|
1523
|
#Input
|
|
2280
|
##Input
|
|
|
1524
|
d = data[:,h]
|
|
2281
|
d = data[:,h]
|
|
|
1525
|
|
|
2282
|
smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
|
|
|
|
|
|
2283
|
signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
|
|
|
|
|
|
2284
|
signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
|
|
|
|
|
|
2285
|
signal0 = signalpn0-n0
|
|
|
|
|
|
2286
|
signal1 = signalpn1-n1
|
|
|
|
|
|
2287
|
snr0 = numpy.sum(signal0/n0)/(nProf-1)
|
|
|
|
|
|
2288
|
snr1 = numpy.sum(signal1/n1)/(nProf-1)
|
|
|
|
|
|
2289
|
#print clean_num_aver[coord,h]
|
|
|
|
|
|
2290
|
if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
|
|
|
1526
|
#Covariance Matrix
|
|
2291
|
#Covariance Matrix
|
|
|
1527
|
D = numpy.diag(d**2/K)
|
|
2292
|
#print h, d.shape
|
|
|
|
|
|
2293
|
D = numpy.diag(d**2)
|
|
|
1528
|
ind = 0
|
|
2294
|
ind = 0
|
|
|
1529
|
for pairs in listComb:
|
|
2295
|
for pairs in listComb:
|
|
|
1530
|
#Coordinates in Covariance Matrix
|
|
2296
|
#Coordinates in Covariance Matrix
|
|
|
@@
-1537,47
+2303,164
class SpectralFitting(Operation):
|
|
|
1537
|
D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
|
|
2303
|
D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
|
|
|
1538
|
D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
|
|
2304
|
D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
|
|
|
1539
|
ind += 1
|
|
2305
|
ind += 1
|
|
|
|
|
|
2306
|
diagD = numpy.zeros(256)
|
|
|
|
|
|
2307
|
if h == 17 :
|
|
|
|
|
|
2308
|
for ii in range(256): diagD[ii] = D[ii,ii]
|
|
|
|
|
|
2309
|
#plt.plot(diagD)
|
|
|
|
|
|
2310
|
#plt.show()
|
|
|
|
|
|
2311
|
|
|
|
|
|
|
2312
|
# print hprint
|
|
|
|
|
|
2313
|
#Dinv=numpy.linalg.inv(D)
|
|
|
|
|
|
2314
|
#L=numpy.linalg.cholesky(Dinv)
|
|
|
|
|
|
2315
|
try:
|
|
|
1540
|
Dinv=numpy.linalg.inv(D)
|
|
2316
|
Dinv=numpy.linalg.inv(D)
|
|
|
1541
|
L=numpy.linalg.cholesky(Dinv)
|
|
2317
|
L=numpy.linalg.cholesky(Dinv)
|
|
|
|
|
|
2318
|
except:
|
|
|
|
|
|
2319
|
Dinv = D*numpy.nan
|
|
|
|
|
|
2320
|
L= D*numpy.nan
|
|
|
1542
|
LT=L.T
|
|
2321
|
LT=L.T
|
|
|
1543
|
|
|
2322
|
|
|
|
1544
|
dp = numpy.dot(LT,d)
|
|
2323
|
dp = numpy.dot(LT,d)
|
|
|
1545
|
|
|
2324
|
|
|
|
1546
|
#Initial values
|
|
2325
|
#Initial values
|
|
|
1547
|
data_spc = self.dataIn.data_spc[coord,:,h]
|
|
2326
|
data_spc = dataOut.data_spc[coord,:,h]
|
|
|
1548
|
|
|
2327
|
|
|
|
1549
|
if (h>0)and(error1[3]<5):
|
|
2328
|
if (h>0)and(error1[3]<5):
|
|
|
1550
|
p0 = self.dataOut.data_param[i,:,h-1]
|
|
2329
|
p0 = dataOut.data_param[i,:,h-1]
|
|
|
1551
|
else:
|
|
2330
|
else:
|
|
|
1552
|
p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
|
|
2331
|
#print("INSIDE ELSE")
|
|
|
1553
|
|
|
2332
|
#print(data_spc.shape,constants,i)
|
|
|
|
|
|
2333
|
p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))# sin el i(data_spc, constants, i)
|
|
|
|
|
|
2334
|
#print ("WAIT_p0",p0)
|
|
|
1554
|
try:
|
|
2335
|
try:
|
|
|
1555
|
#Least Squares
|
|
2336
|
#Least Squares
|
|
|
|
|
|
2337
|
#print (dp,LT,constants)
|
|
|
|
|
|
2338
|
#value =self.__residFunction(p0,dp,LT,constants)
|
|
|
|
|
|
2339
|
#print ("valueREADY",value.shape, type(value))
|
|
|
|
|
|
2340
|
#optimize.leastsq(value)
|
|
|
1556
|
minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
|
|
2341
|
minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
|
|
|
|
|
|
2342
|
|
|
|
|
|
|
2343
|
# print(minp)
|
|
|
1557
|
#minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
|
|
2344
|
#minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
|
|
|
1558
|
#Chi square error
|
|
2345
|
#Chi square error
|
|
|
|
|
|
2346
|
#print(minp,covp.infodict,mesg,ier)
|
|
|
|
|
|
2347
|
#print("REALIZA OPTIMIZ")
|
|
|
1559
|
error0 = numpy.sum(infodict['fvec']**2)/(2*N)
|
|
2348
|
error0 = numpy.sum(infodict['fvec']**2)/(2*N)
|
|
|
1560
|
#Error with Jacobian
|
|
2349
|
#Error with Jacobian
|
|
|
1561
|
error1 = self.dataOut.library.errorFunction(minp,constants,LT)
|
|
2350
|
error1 = self.library.errorFunction(minp,constants,LT)
|
|
|
|
|
|
2351
|
# print self.__residFunction(p0,dp,LT, constants)
|
|
|
|
|
|
2352
|
# print infodict['fvec']
|
|
|
|
|
|
2353
|
# print self.__residFunction(minp,dp,LT,constants)
|
|
|
|
|
|
2354
|
|
|
|
1562
|
except:
|
|
2355
|
except:
|
|
|
1563
|
minp = p0*numpy.nan
|
|
2356
|
minp = p0*numpy.nan
|
|
|
1564
|
error0 = numpy.nan
|
|
2357
|
error0 = numpy.nan
|
|
|
1565
|
error1 = p0*numpy.nan
|
|
2358
|
error1 = p0*numpy.nan
|
|
|
1566
|
|
|
2359
|
#print ("EXCEPT 0000000000")
|
|
|
|
|
|
2360
|
# s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0))
|
|
|
|
|
|
2361
|
# covp = covp*s_sq
|
|
|
|
|
|
2362
|
# #print("TRY___________________________________________1")
|
|
|
|
|
|
2363
|
# error = []
|
|
|
|
|
|
2364
|
# for ip in range(len(minp)):
|
|
|
|
|
|
2365
|
# try:
|
|
|
|
|
|
2366
|
# error.append(numpy.absolute(covp[ip][ip])**0.5)
|
|
|
|
|
|
2367
|
# except:
|
|
|
|
|
|
2368
|
# error.append( 0.00 )
|
|
|
|
|
|
2369
|
else :
|
|
|
|
|
|
2370
|
data_spc = dataOut.data_spc[coord,:,h]
|
|
|
|
|
|
2371
|
p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
|
|
|
|
|
|
2372
|
minp = p0*numpy.nan
|
|
|
|
|
|
2373
|
error0 = numpy.nan
|
|
|
|
|
|
2374
|
error1 = p0*numpy.nan
|
|
|
1567
|
#Save
|
|
2375
|
#Save
|
|
|
1568
|
if self.dataOut.data_param is None:
|
|
2376
|
if dataOut.data_param is None:
|
|
|
1569
|
self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
|
|
2377
|
dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
|
|
|
1570
|
self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
|
|
2378
|
dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
|
|
|
|
|
|
2379
|
|
|
|
|
|
|
2380
|
dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
|
|
|
|
|
|
2381
|
dataOut.data_param[i,:,h] = minp
|
|
|
|
|
|
2382
|
#print(minp)
|
|
|
|
|
|
2383
|
#print("FIN")
|
|
|
|
|
|
2384
|
#print ("DATA",minp.shape)
|
|
|
|
|
|
2385
|
|
|
|
|
|
|
2386
|
#plt.plot(dataOut.data_param[0,3,:])
|
|
|
|
|
|
2387
|
#print(dataOut.data_param[:,3,:])
|
|
|
|
|
|
2388
|
#dataOut.data_errorC = numpy.zeros((nGroups, nHeights,1))*numpy.nan
|
|
|
|
|
|
2389
|
for ht in range(nHeights-1) :
|
|
|
|
|
|
2390
|
smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
|
|
|
|
|
|
2391
|
dataOut.data_paramC[4*i,ht,1] = smooth
|
|
|
|
|
|
2392
|
signalpn0 = (coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
|
|
|
|
|
|
2393
|
signalpn1 = (coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
|
|
|
|
|
|
2394
|
|
|
|
|
|
|
2395
|
#val0 = WHERE(signalpn0 > 0,cval0)
|
|
|
|
|
|
2396
|
val0 = (signalpn0 > 0).nonzero()
|
|
|
|
|
|
2397
|
val0 = val0[0]
|
|
|
|
|
|
2398
|
#print('hvalid:',hvalid)
|
|
|
|
|
|
2399
|
#print('valid', valid)
|
|
|
|
|
|
2400
|
if len(val0) == 0 : val0_npoints = nProf
|
|
|
|
|
|
2401
|
else : val0_npoints = len(val0)
|
|
|
|
|
|
2402
|
|
|
|
|
|
|
2403
|
#val1 = WHERE(signalpn1 > 0,cval1)
|
|
|
|
|
|
2404
|
val1 = (signalpn1 > 0).nonzero()
|
|
|
|
|
|
2405
|
val1 = val1[0]
|
|
|
|
|
|
2406
|
if len(val1) == 0 : val1_npoints = nProf
|
|
|
|
|
|
2407
|
else : val1_npoints = len(val1)
|
|
|
|
|
|
2408
|
|
|
|
|
|
|
2409
|
dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
|
|
|
|
|
|
2410
|
dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
|
|
|
|
|
|
2411
|
|
|
|
|
|
|
2412
|
signal0 = (signalpn0-n0) # > 0
|
|
|
|
|
|
2413
|
vali = (signal0 < 0).nonzero()
|
|
|
|
|
|
2414
|
vali = vali[0]
|
|
|
|
|
|
2415
|
if len(vali) > 0 : signal0[vali] = 0
|
|
|
|
|
|
2416
|
signal1 = (signalpn1-n1) #> 0
|
|
|
|
|
|
2417
|
vali = (signal1 < 0).nonzero()
|
|
|
|
|
|
2418
|
vali = vali[0]
|
|
|
|
|
|
2419
|
if len(vali) > 0 : signal1[vali] = 0
|
|
|
|
|
|
2420
|
snr0 = numpy.sum(signal0/n0)/(nProf-1)
|
|
|
|
|
|
2421
|
snr1 = numpy.sum(signal1/n1)/(nProf-1)
|
|
|
|
|
|
2422
|
doppler = absc[1:]
|
|
|
|
|
|
2423
|
if snr0 >= snrth and snr1 >= snrth and smooth :
|
|
|
|
|
|
2424
|
signalpn0_n0 = signalpn0
|
|
|
|
|
|
2425
|
signalpn0_n0[val0] = signalpn0[val0] - n0
|
|
|
|
|
|
2426
|
mom0 = self.moments(doppler,signalpn0-n0,nProf)
|
|
|
|
|
|
2427
|
# sigtmp= numpy.transpose(numpy.tile(signalpn0, [4,1]))
|
|
|
|
|
|
2428
|
# momt= self.__calculateMoments( sigtmp, doppler , n0 )
|
|
|
|
|
|
2429
|
signalpn1_n1 = signalpn1
|
|
|
|
|
|
2430
|
signalpn1_n1[val1] = signalpn1[val1] - n1
|
|
|
|
|
|
2431
|
mom1 = self.moments(doppler,signalpn1_n1,nProf)
|
|
|
|
|
|
2432
|
dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
|
|
|
|
|
|
2433
|
dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
|
|
|
|
|
|
2434
|
# if graph == 1 :
|
|
|
|
|
|
2435
|
# window, 13
|
|
|
|
|
|
2436
|
# plot,doppler,signalpn0
|
|
|
|
|
|
2437
|
# oplot,doppler,signalpn1,linest=1
|
|
|
|
|
|
2438
|
# oplot,mom0(0)*doppler/doppler,signalpn0
|
|
|
|
|
|
2439
|
# oplot,mom1(0)*doppler/doppler,signalpn1
|
|
|
|
|
|
2440
|
# print,interval/12.,beam,45+ht*15,snr0,snr1,mom0(0),mom1(0),mom0(1),mom1(1)
|
|
|
|
|
|
2441
|
#ENDIF
|
|
|
|
|
|
2442
|
#ENDIF
|
|
|
|
|
|
2443
|
#ENDFOR End height
|
|
|
|
|
|
2444
|
#plt.show()
|
|
|
|
|
|
2445
|
#print dataOut.data_param[i,3,:]
|
|
|
|
|
|
2446
|
# if self.__dataReady:
|
|
|
|
|
|
2447
|
# dataOut.flagNoData = False
|
|
|
|
|
|
2448
|
#print dataOut.data_error[:,3,:]
|
|
|
|
|
|
2449
|
dataOut.data_spc = jspectra
|
|
|
|
|
|
2450
|
if getSNR:
|
|
|
|
|
|
2451
|
listChannels = groupArray.reshape((groupArray.size))
|
|
|
|
|
|
2452
|
listChannels.sort()
|
|
|
|
|
|
2453
|
|
|
|
|
|
|
2454
|
dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels])
|
|
|
|
|
|
2455
|
return dataOut
|
|
|
1571
|
|
|
2456
|
|
|
|
1572
|
self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
|
|
|
|
|
|
1573
|
self.dataOut.data_param[i,:,h] = minp
|
|
|
|
|
|
1574
|
return
|
|
|
|
|
|
1575
|
|
|
2457
|
|
|
|
1576
|
def __residFunction(self, p, dp, LT, constants):
|
|
2458
|
def __residFunction(self, p, dp, LT, constants):
|
|
|
1577
|
|
|
2459
|
|
|
|
1578
|
fm = self.dataOut.library.modelFunction(p, constants)
|
|
2460
|
fm = self.library.modelFunction(p, constants)
|
|
|
1579
|
fmp=numpy.dot(LT,fm)
|
|
2461
|
fmp=numpy.dot(LT,fm)
|
|
|
1580
|
|
|
2462
|
#print ("DONE",dp -fmp)
|
|
|
|
|
|
2463
|
#print ("ok")
|
|
|
1581
|
return dp-fmp
|
|
2464
|
return dp-fmp
|
|
|
1582
|
|
|
2465
|
|
|
|
1583
|
def __getSNR(self, z, noise):
|
|
2466
|
def __getSNR(self, z, noise):
|
|
|
@@
-1587,7
+2470,7
class SpectralFitting(Operation):
|
|
|
1587
|
SNR = SNR.T
|
|
2470
|
SNR = SNR.T
|
|
|
1588
|
return SNR
|
|
2471
|
return SNR
|
|
|
1589
|
|
|
2472
|
|
|
|
1590
|
def __chisq(p,chindex,hindex):
|
|
2473
|
def __chisq(self,p,chindex,hindex):
|
|
|
1591
|
#similar to Resid but calculates CHI**2
|
|
2474
|
#similar to Resid but calculates CHI**2
|
|
|
1592
|
[LT,d,fm]=setupLTdfm(p,chindex,hindex)
|
|
2475
|
[LT,d,fm]=setupLTdfm(p,chindex,hindex)
|
|
|
1593
|
dp=numpy.dot(LT,d)
|
|
2476
|
dp=numpy.dot(LT,d)
|
|
|
@@
-2152,7
+3035,7
class WindProfiler(Operation):
|
|
|
2152
|
absc = dataOut.abscissaList[:-1]
|
|
3035
|
absc = dataOut.abscissaList[:-1]
|
|
|
2153
|
# noise = dataOut.noise
|
|
3036
|
# noise = dataOut.noise
|
|
|
2154
|
heightList = dataOut.heightList
|
|
3037
|
heightList = dataOut.heightList
|
|
|
2155
|
SNR = dataOut.data_snr
|
|
3038
|
SNR = dataOut.data_SNR
|
|
|
2156
|
|
|
3039
|
|
|
|
2157
|
if technique == 'DBS':
|
|
3040
|
if technique == 'DBS':
|
|
|
2158
|
|
|
3041
|
|
|
|
@@
-2160,7
+3043,7
class WindProfiler(Operation):
|
|
|
2160
|
kwargs['heightList'] = heightList
|
|
3043
|
kwargs['heightList'] = heightList
|
|
|
2161
|
kwargs['SNR'] = SNR
|
|
3044
|
kwargs['SNR'] = SNR
|
|
|
2162
|
|
|
3045
|
|
|
|
2163
|
dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
|
|
3046
|
dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
|
|
|
2164
|
dataOut.utctimeInit = dataOut.utctime
|
|
3047
|
dataOut.utctimeInit = dataOut.utctime
|
|
|
2165
|
dataOut.outputInterval = dataOut.paramInterval
|
|
3048
|
dataOut.outputInterval = dataOut.paramInterval
|
|
|
2166
|
|
|
3049
|
|
|
|
@@
-2337,13
+3220,23
class EWDriftsEstimation(Operation):
|
|
|
2337
|
for i in rango:
|
|
3220
|
for i in rango:
|
|
|
2338
|
x = heiRang*math.cos(phi[i])
|
|
3221
|
x = heiRang*math.cos(phi[i])
|
|
|
2339
|
y1 = velRadial[i,:]
|
|
3222
|
y1 = velRadial[i,:]
|
|
|
2340
|
f1 = interpolate.interp1d(x,y1,kind = 'cubic')
|
|
3223
|
vali= (numpy.isfinite(y1)==True).nonzero()
|
|
|
|
|
|
3224
|
y1=y1[vali]
|
|
|
|
|
|
3225
|
x = x[vali]
|
|
|
|
|
|
3226
|
f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
|
|
|
2341
|
|
|
3227
|
|
|
|
|
|
|
3228
|
#heiRang1 = x*math.cos(phi[maxid])
|
|
|
2342
|
x1 = heiRang1
|
|
3229
|
x1 = heiRang1
|
|
|
2343
|
y11 = f1(x1)
|
|
3230
|
y11 = f1(x1)
|
|
|
2344
|
|
|
3231
|
|
|
|
2345
|
y2 = SNR[i,:]
|
|
3232
|
y2 = SNR[i,:]
|
|
|
2346
|
f2 = interpolate.interp1d(x,y2,kind = 'cubic')
|
|
3233
|
#print 'snr ', y2
|
|
|
|
|
|
3234
|
x = heiRang*math.cos(phi[i])
|
|
|
|
|
|
3235
|
vali= (y2 != -1).nonzero()
|
|
|
|
|
|
3236
|
y2 = y2[vali]
|
|
|
|
|
|
3237
|
x = x[vali]
|
|
|
|
|
|
3238
|
#print 'snr ',y2
|
|
|
|
|
|
3239
|
f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
|
|
|
2347
|
y21 = f2(x1)
|
|
3240
|
y21 = f2(x1)
|
|
|
2348
|
|
|
3241
|
|
|
|
2349
|
velRadial1[i,:] = y11
|
|
3242
|
velRadial1[i,:] = y11
|
|
|
@@
-2351,35
+3244,185
class EWDriftsEstimation(Operation):
|
|
|
2351
|
|
|
3244
|
|
|
|
2352
|
return heiRang1, velRadial1, SNR1
|
|
3245
|
return heiRang1, velRadial1, SNR1
|
|
|
2353
|
|
|
3246
|
|
|
|
|
|
|
3247
|
|
|
|
|
|
|
3248
|
|
|
|
2354
|
def run(self, dataOut, zenith, zenithCorrection):
|
|
3249
|
def run(self, dataOut, zenith, zenithCorrection):
|
|
|
|
|
|
3250
|
import matplotlib.pyplot as plt
|
|
|
2355
|
heiRang = dataOut.heightList
|
|
3251
|
heiRang = dataOut.heightList
|
|
|
2356
|
velRadial = dataOut.data_param[:,3,:]
|
|
3252
|
velRadial = dataOut.data_param[:,3,:]
|
|
|
2357
|
SNR = dataOut.data_snr
|
|
3253
|
velRadialm = dataOut.data_param[:,2:4,:]*-1
|
|
|
2358
|
|
|
3254
|
|
|
|
|
|
|
3255
|
rbufc=dataOut.data_paramC[:,:,0]
|
|
|
|
|
|
3256
|
ebufc=dataOut.data_paramC[:,:,1]
|
|
|
|
|
|
3257
|
SNR = dataOut.data_SNR
|
|
|
|
|
|
3258
|
velRerr = dataOut.data_error[:,4,:]
|
|
|
|
|
|
3259
|
moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
|
|
|
|
|
|
3260
|
dataOut.moments=moments
|
|
|
|
|
|
3261
|
# Coherent
|
|
|
|
|
|
3262
|
smooth_wC = ebufc[0,:]
|
|
|
|
|
|
3263
|
p_w0C = rbufc[0,:]
|
|
|
|
|
|
3264
|
p_w1C = rbufc[1,:]
|
|
|
|
|
|
3265
|
w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
|
|
|
|
|
|
3266
|
t_wC = rbufc[3,:]
|
|
|
|
|
|
3267
|
my_nbeams = 2
|
|
|
|
|
|
3268
|
|
|
|
|
|
|
3269
|
# plt.plot(w_wC)
|
|
|
|
|
|
3270
|
# plt.show()
|
|
|
2359
|
zenith = numpy.array(zenith)
|
|
3271
|
zenith = numpy.array(zenith)
|
|
|
2360
|
zenith -= zenithCorrection
|
|
3272
|
zenith -= zenithCorrection
|
|
|
2361
|
zenith *= numpy.pi/180
|
|
3273
|
zenith *= numpy.pi/180
|
|
|
2362
|
|
|
3274
|
if zenithCorrection != 0 :
|
|
|
2363
|
heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
|
|
3275
|
heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
|
|
|
|
|
|
3276
|
else :
|
|
|
|
|
|
3277
|
heiRang1 = heiRang
|
|
|
|
|
|
3278
|
velRadial1 = velRadial
|
|
|
|
|
|
3279
|
SNR1 = SNR
|
|
|
2364
|
|
|
3280
|
|
|
|
2365
|
alp = zenith[0]
|
|
3281
|
alp = zenith[0]
|
|
|
2366
|
bet = zenith[1]
|
|
3282
|
bet = zenith[1]
|
|
|
2367
|
|
|
3283
|
|
|
|
|
|
|
3284
|
#t_w(bad) = t_wC(bad)
|
|
|
|
|
|
3285
|
#t_w_err(bad)=!values.f_nan
|
|
|
|
|
|
3286
|
|
|
|
2368
|
w_w = velRadial1[0,:]
|
|
3287
|
w_w = velRadial1[0,:]
|
|
|
2369
|
w_e = velRadial1[1,:]
|
|
3288
|
w_e = velRadial1[1,:]
|
|
|
|
|
|
3289
|
w_w_err = velRerr[0,:]
|
|
|
|
|
|
3290
|
w_e_err = velRerr[1,:]
|
|
|
|
|
|
3291
|
#plt.plot(w_w)
|
|
|
|
|
|
3292
|
#plt.show()
|
|
|
|
|
|
3293
|
#plt.plot(w_e)
|
|
|
|
|
|
3294
|
#plt.show()
|
|
|
|
|
|
3295
|
# bad = where((chisq_w GT 2.5 AND abs(w_w_err) GT 1. AND finite(w_wC))
|
|
|
|
|
|
3296
|
# OR abs(w_w) GT 200. OR (NOT finite(w_w)-254) OR ABS(w_w_err) GT 100, cbad)
|
|
|
|
|
|
3297
|
val = (numpy.isfinite(w_w)==False).nonzero()
|
|
|
|
|
|
3298
|
val = val[0]
|
|
|
|
|
|
3299
|
bad = val
|
|
|
|
|
|
3300
|
if len(bad) > 0 :
|
|
|
|
|
|
3301
|
w_w[bad] = w_wC[bad]
|
|
|
|
|
|
3302
|
w_w_err[bad]= numpy.nan
|
|
|
|
|
|
3303
|
if my_nbeams == 2:
|
|
|
|
|
|
3304
|
smooth_eC=ebufc[4,:]
|
|
|
|
|
|
3305
|
p_e0C = rbufc[4,:]
|
|
|
|
|
|
3306
|
p_e1C = rbufc[5,:]
|
|
|
|
|
|
3307
|
w_eC = rbufc[6,:]*-1
|
|
|
|
|
|
3308
|
t_eC = rbufc[7,:]
|
|
|
|
|
|
3309
|
val = (numpy.isfinite(w_e)==False).nonzero()
|
|
|
|
|
|
3310
|
val = val[0]
|
|
|
|
|
|
3311
|
bad = val
|
|
|
|
|
|
3312
|
if len(bad) > 0 :
|
|
|
|
|
|
3313
|
w_e[bad] = w_eC[bad]
|
|
|
|
|
|
3314
|
w_e_err[bad]= numpy.nan
|
|
|
2370
|
|
|
3315
|
|
|
|
2371
|
w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
|
|
3316
|
w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
|
|
|
2372
|
u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
|
|
3317
|
u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
|
|
|
|
|
|
3318
|
#plt.plot(w)
|
|
|
|
|
|
3319
|
#plt.show()
|
|
|
|
|
|
3320
|
#error
|
|
|
|
|
|
3321
|
w_err = numpy.sqrt((w_w_err*numpy.sin(bet))**2.+(w_e_err*numpy.sin(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
|
|
|
|
|
|
3322
|
u_err = numpy.sqrt((w_w_err*numpy.cos(bet))**2.+(w_e_err*numpy.cos(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
|
|
|
2373
|
|
|
3323
|
|
|
|
2374
|
winds = numpy.vstack((u,w))
|
|
3324
|
winds = numpy.vstack((w,u))
|
|
|
2375
|
|
|
3325
|
|
|
|
2376
|
dataOut.heightList = heiRang1
|
|
3326
|
dataOut.heightList = heiRang1
|
|
|
2377
|
dataOut.data_output = winds
|
|
3327
|
dataOut.data_output = winds
|
|
|
2378
|
dataOut.data_snr = SNR1
|
|
3328
|
#dataOut.data_SNR = SNR1
|
|
|
2379
|
|
|
3329
|
snr1 = 10*numpy.log10(SNR1[0])
|
|
|
|
|
|
3330
|
dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
|
|
|
2380
|
dataOut.utctimeInit = dataOut.utctime
|
|
3331
|
dataOut.utctimeInit = dataOut.utctime
|
|
|
2381
|
dataOut.outputInterval = dataOut.timeInterval
|
|
3332
|
dataOut.outputInterval = dataOut.timeInterval
|
|
|
2382
|
return
|
|
3333
|
|
|
|
|
|
|
3334
|
hei_aver0 = 218
|
|
|
|
|
|
3335
|
jrange = 450 #900 para HA drifts
|
|
|
|
|
|
3336
|
deltah = 15.0 #dataOut.spacing(0)
|
|
|
|
|
|
3337
|
h0 = 0.0 #dataOut.first_height(0)
|
|
|
|
|
|
3338
|
heights = dataOut.heightList
|
|
|
|
|
|
3339
|
nhei = len(heights)
|
|
|
|
|
|
3340
|
|
|
|
|
|
|
3341
|
range1 = numpy.arange(nhei) * deltah + h0
|
|
|
|
|
|
3342
|
|
|
|
|
|
|
3343
|
#jhei = WHERE(range1 GE hei_aver0 , jcount)
|
|
|
|
|
|
3344
|
jhei = (range1 >= hei_aver0).nonzero()
|
|
|
|
|
|
3345
|
if len(jhei[0]) > 0 :
|
|
|
|
|
|
3346
|
h0_index = jhei[0][0] # Initial height for getting averages 218km
|
|
|
|
|
|
3347
|
|
|
|
|
|
|
3348
|
mynhei = 7
|
|
|
|
|
|
3349
|
nhei_avg = int(jrange/deltah)
|
|
|
|
|
|
3350
|
h_avgs = int(nhei_avg/mynhei)
|
|
|
|
|
|
3351
|
nhei_avg = h_avgs*(mynhei-1)+mynhei
|
|
|
|
|
|
3352
|
|
|
|
|
|
|
3353
|
navgs = numpy.zeros(mynhei,dtype='float')
|
|
|
|
|
|
3354
|
delta_h = numpy.zeros(mynhei,dtype='float')
|
|
|
|
|
|
3355
|
range_aver = numpy.zeros(mynhei,dtype='float')
|
|
|
|
|
|
3356
|
for ih in range( mynhei-1 ):
|
|
|
|
|
|
3357
|
range_aver[ih] = numpy.sum(range1[h0_index+h_avgs*ih:h0_index+h_avgs*(ih+1)-0])/h_avgs
|
|
|
|
|
|
3358
|
navgs[ih] = h_avgs
|
|
|
|
|
|
3359
|
delta_h[ih] = deltah*h_avgs
|
|
|
|
|
|
3360
|
|
|
|
|
|
|
3361
|
range_aver[mynhei-1] = numpy.sum(range1[h0_index:h0_index+6*h_avgs-0])/(6*h_avgs)
|
|
|
|
|
|
3362
|
navgs[mynhei-1] = 6*h_avgs
|
|
|
|
|
|
3363
|
delta_h[mynhei-1] = deltah*6*h_avgs
|
|
|
|
|
|
3364
|
|
|
|
|
|
|
3365
|
wA = w[h0_index:h0_index+nhei_avg-0]
|
|
|
|
|
|
3366
|
wA_err = w_err[h0_index:h0_index+nhei_avg-0]
|
|
|
|
|
|
3367
|
#print(wA, wA_err)
|
|
|
|
|
|
3368
|
for i in range(5) :
|
|
|
|
|
|
3369
|
vals = wA[i*h_avgs:(i+1)*h_avgs-0]
|
|
|
|
|
|
3370
|
errs = wA_err[i*h_avgs:(i+1)*h_avgs-0]
|
|
|
|
|
|
3371
|
avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2.)
|
|
|
|
|
|
3372
|
sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
|
|
|
|
|
|
3373
|
wA[6*h_avgs+i] = avg
|
|
|
|
|
|
3374
|
wA_err[6*h_avgs+i] = sigma
|
|
|
|
|
|
3375
|
|
|
|
|
|
|
3376
|
|
|
|
|
|
|
3377
|
vals = wA[0:6*h_avgs-0]
|
|
|
|
|
|
3378
|
errs=wA_err[0:6*h_avgs-0]
|
|
|
|
|
|
3379
|
avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2)
|
|
|
|
|
|
3380
|
sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
|
|
|
|
|
|
3381
|
wA[nhei_avg-1] = avg
|
|
|
|
|
|
3382
|
wA_err[nhei_avg-1] = sigma
|
|
|
|
|
|
3383
|
|
|
|
|
|
|
3384
|
wA = wA[6*h_avgs:nhei_avg-0]
|
|
|
|
|
|
3385
|
wA_err=wA_err[6*h_avgs:nhei_avg-0]
|
|
|
|
|
|
3386
|
if my_nbeams == 2 :
|
|
|
|
|
|
3387
|
|
|
|
|
|
|
3388
|
uA = u[h0_index:h0_index+nhei_avg]
|
|
|
|
|
|
3389
|
uA_err=u_err[h0_index:h0_index+nhei_avg]
|
|
|
|
|
|
3390
|
|
|
|
|
|
|
3391
|
for i in range(5) :
|
|
|
|
|
|
3392
|
vals = uA[i*h_avgs:(i+1)*h_avgs-0]
|
|
|
|
|
|
3393
|
errs=uA_err[i*h_avgs:(i+1)*h_avgs-0]
|
|
|
|
|
|
3394
|
avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2.)
|
|
|
|
|
|
3395
|
sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
|
|
|
|
|
|
3396
|
uA[6*h_avgs+i] = avg
|
|
|
|
|
|
3397
|
uA_err[6*h_avgs+i]=sigma
|
|
|
|
|
|
3398
|
|
|
|
|
|
|
3399
|
vals = uA[0:6*h_avgs-0]
|
|
|
|
|
|
3400
|
errs = uA_err[0:6*h_avgs-0]
|
|
|
|
|
|
3401
|
avg = numpy.nansum(vals/errs**2.)/numpy.nansum(1./errs**2.)
|
|
|
|
|
|
3402
|
sigma = numpy.sqrt(1./numpy.nansum(1./errs**2.))
|
|
|
|
|
|
3403
|
uA[nhei_avg-1] = avg
|
|
|
|
|
|
3404
|
uA_err[nhei_avg-1] = sigma
|
|
|
|
|
|
3405
|
uA = uA[6*h_avgs:nhei_avg-0]
|
|
|
|
|
|
3406
|
uA_err = uA_err[6*h_avgs:nhei_avg-0]
|
|
|
|
|
|
3407
|
|
|
|
|
|
|
3408
|
dataOut.drifts_avg = numpy.vstack((wA,uA))
|
|
|
|
|
|
3409
|
#print(dataOut.drifts_avg)
|
|
|
|
|
|
3410
|
tini=time.localtime(dataOut.utctime)
|
|
|
|
|
|
3411
|
datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
|
|
|
|
|
|
3412
|
nfile = '/home/pcondor/Database/ewdriftsschain2019/jro'+datefile+'drifts_sch3.txt'
|
|
|
|
|
|
3413
|
#print(dataOut.drifts_avg)
|
|
|
|
|
|
3414
|
f1 = open(nfile,'a')
|
|
|
|
|
|
3415
|
#print(nfile)
|
|
|
|
|
|
3416
|
#f.write(datefile)
|
|
|
|
|
|
3417
|
#numpy.savetxt(f,[datefile,datefile],fmt='%10s')
|
|
|
|
|
|
3418
|
datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
|
|
|
|
|
|
3419
|
driftavgstr=str(dataOut.drifts_avg)
|
|
|
|
|
|
3420
|
#f1.write(datedriftavg)
|
|
|
|
|
|
3421
|
#f1.write(driftavgstr)
|
|
|
|
|
|
3422
|
numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
|
|
|
|
|
|
3423
|
numpy.savetxt(f1,dataOut.drifts_avg,fmt='%10.2f')
|
|
|
|
|
|
3424
|
f1.close()
|
|
|
|
|
|
3425
|
return dataOut
|
|
|
2383
|
|
|
3426
|
|
|
|
2384
|
#--------------- Non Specular Meteor ----------------
|
|
3427
|
#--------------- Non Specular Meteor ----------------
|
|
|
2385
|
|
|
3428
|
|
|
|
@@
-2656,7
+3699,7
class SMDetection(Operation):
|
|
|
2656
|
channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
|
|
3699
|
channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
|
|
|
2657
|
meteorOps = SMOperations()
|
|
3700
|
meteorOps = SMOperations()
|
|
|
2658
|
pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
|
|
3701
|
pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
|
|
|
2659
|
heiRang = dataOut.heightList
|
|
3702
|
heiRang = dataOut.getHeiRange()
|
|
|
2660
|
#Get Beacon signal - No Beacon signal anymore
|
|
3703
|
#Get Beacon signal - No Beacon signal anymore
|
|
|
2661
|
# newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
|
|
3704
|
# newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
|
|
|
2662
|
#
|
|
3705
|
#
|
|
|
@@
-2718,7
+3761,7
class SMDetection(Operation):
|
|
|
2718
|
|
|
3761
|
|
|
|
2719
|
#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
|
|
3762
|
#************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
|
|
|
2720
|
#Parameters
|
|
3763
|
#Parameters
|
|
|
2721
|
heiRange = dataOut.heightList
|
|
3764
|
heiRange = dataOut.getHeiRange()
|
|
|
2722
|
rangeInterval = heiRange[1] - heiRange[0]
|
|
3765
|
rangeInterval = heiRange[1] - heiRange[0]
|
|
|
2723
|
rangeLimit = multDet_rangeLimit/rangeInterval
|
|
3766
|
rangeLimit = multDet_rangeLimit/rangeInterval
|
|
|
2724
|
timeLimit = multDet_timeLimit/dataOut.timeInterval
|
|
3767
|
timeLimit = multDet_timeLimit/dataOut.timeInterval
|
|
|
@@
-3884,3
+4927,4
class SMOperations():
|
|
|
3884
|
# error[indInvalid1] = 13
|
|
4927
|
# error[indInvalid1] = 13
|
|
|
3885
|
#
|
|
4928
|
#
|
|
|
3886
|
# return heights, error
|
|
4929
|
# return heights, error
|
|
|
|
|
|
4930
|
|
