##// END OF EJS Templates
schain
RSS
Drifts tested
joabAM -
r1740:7f5b085e2124
parent child
Show More
Side by Side Unified
Context file:
r1740:7f5b085e2124 Loading diff...:
Show file before | Show file after | Hide comments
@@ -459,6 +459,9 class Voltage(JROData):
459 powerdB = numpy.squeeze(powerdB)
459 powerdB = numpy.squeeze(powerdB)
460
460
461 return powerdB
461 return powerdB
462 @property
463 def data_pow(self):
464 return self.getPower()
462
465
463 @property
466 @property
464 def timeInterval(self):
467 def timeInterval(self):
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -65,9 +65,13 class SpectraPlot(Plot):
65 self.update_list(dataOut)
65 self.update_list(dataOut)
66 data = {}
66 data = {}
67 meta = {}
67 meta = {}
68
69 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
69 if dataOut.type == "Parameters":
70 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
71 spc = 10*numpy.log10(dataOut.data_spc/(dataOut.nProfiles))
72 else:
70 noise = 10*numpy.log10(dataOut.getNoise()/norm)
73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74
71 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
72 for ch in range(dataOut.nChannels):
76 for ch in range(dataOut.nChannels):
73 if hasattr(dataOut.normFactor,'ndim'):
77 if hasattr(dataOut.normFactor,'ndim'):
@@ -725,7 +729,7 class RTIPlot(Plot):
725 )
729 )
726 if self.showprofile:
730 if self.showprofile:
727 ax.plot_profile = self.pf_axes[n].plot(
731 ax.plot_profile = self.pf_axes[n].plot(
728 data['rti'][n], self.y)[0]
732 data[self.CODE][n], self.y)[0]
729 if "noise" in self.data:
733 if "noise" in self.data:
730 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
731 color="k", linestyle="dashed", lw=1)[0]
735 color="k", linestyle="dashed", lw=1)[0]
@@ -737,7 +741,7 class RTIPlot(Plot):
737 cmap=plt.get_cmap(self.colormap)
741 cmap=plt.get_cmap(self.colormap)
738 )
742 )
739 if self.showprofile:
743 if self.showprofile:
740 ax.plot_profile.set_data(data['rti'][n], self.y)
744 ax.plot_profile.set_data(data[self.CODE][n], self.y)
741 if "noise" in self.data:
745 if "noise" in self.data:
742 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
743 color="k", linestyle="dashed", lw=1)[0]
747 color="k", linestyle="dashed", lw=1)[0]
@@ -893,12 +897,15 class NoisePlot(Plot):
893 self.titles = ['Noise']
897 self.titles = ['Noise']
894 self.colorbar = False
898 self.colorbar = False
895 self.plots_adjust.update({'right': 0.85 })
899 self.plots_adjust.update({'right': 0.85 })
900 self.titles = ['Noise Plot']
896
901
897 def update(self, dataOut):
902 def update(self, dataOut):
898
903
899 data = {}
904 data = {}
900 meta = {}
905 meta = {}
901 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
906 noise = 10*numpy.log10(dataOut.getNoise())
907 noise = noise.reshape(dataOut.nChannels, 1)
908 data['noise'] = noise
902 meta['yrange'] = numpy.array([])
909 meta['yrange'] = numpy.array([])
903
910
904 return data, meta
911 return data, meta
Show file before | Show file after | Hide comments
@@ -445,6 +445,7 class HDFWriter(Operation):
445 path = None
445 path = None
446 setFile = None
446 setFile = None
447 fp = None
447 fp = None
448 ds = None
448 firsttime = True
449 firsttime = True
449 #Configurations
450 #Configurations
450 blocksPerFile = None
451 blocksPerFile = None
@@ -490,6 +491,8 class HDFWriter(Operation):
490 if self.metadataList is None:
491 if self.metadataList is None:
491 self.metadataList = self.dataOut.metadata_list
492 self.metadataList = self.dataOut.metadata_list
492
493
494 self.metadataList = list(set(self.metadataList))
495
493 tableList = []
496 tableList = []
494 dsList = []
497 dsList = []
495
498
@@ -504,7 +507,7 class HDFWriter(Operation):
504
507
505 if dataAux is None:
508 if dataAux is None:
506 continue
509 continue
507 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
510 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float_)):
508 dsDict['nDim'] = 0
511 dsDict['nDim'] = 0
509 else:
512 else:
510 dsDict['nDim'] = len(dataAux.shape)
513 dsDict['nDim'] = len(dataAux.shape)
@@ -514,6 +517,7 class HDFWriter(Operation):
514
517
515 dsList.append(dsDict)
518 dsList.append(dsDict)
516
519
520 self.blockIndex = 0
517 self.dsList = dsList
521 self.dsList = dsList
518 self.currentDay = self.dataOut.datatime.date()
522 self.currentDay = self.dataOut.datatime.date()
519
523
@@ -763,7 +767,7 class HDFWriter(Operation):
763 self.ds = dtsets
767 self.ds = dtsets
764 self.data = data
768 self.data = data
765 self.firsttime = True
769 self.firsttime = True
766 self.blockIndex = 0
770
767 return
771 return
768
772
769 def putData(self):
773 def putData(self):
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (698 lines changed) Show them Hide them
@@ -24,7 +24,7 import warnings
24 from numpy import NaN
24 from numpy import NaN
25 from scipy.optimize.optimize import OptimizeWarning
25 from scipy.optimize.optimize import OptimizeWarning
26 warnings.filterwarnings('ignore')
26 warnings.filterwarnings('ignore')
27
27 import json
28 import os
28 import os
29 import csv
29 import csv
30 from scipy import signal
30 from scipy import signal
@@ -3176,7 +3176,6 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
3176 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
3176 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
3177 return y
3177 return y
3178
3178
3179
3180 class SpectralFitting(Operation):
3179 class SpectralFitting(Operation):
3181 '''
3180 '''
3182 Function GetMoments()
3181 Function GetMoments()
@@ -3189,72 +3188,16 class SpectralFitting(Operation):
3189 __dataReady = False
3188 __dataReady = False
3190 bloques = None
3189 bloques = None
3191 bloque0 = None
3190 bloque0 = None
3192 index = 0
3193 fint = 0
3194 buffer = 0
3195 buffer2 = 0
3196 buffer3 = 0
3197
3191
3198 def __init__(self):
3192 def __init__(self):
3199 Operation.__init__(self)
3193 Operation.__init__(self)
3200 self.i=0
3194 self.i=0
3201 self.isConfig = False
3195 self.isConfig = False
3202
3196
3203
3197 def setup(self,nChan,nProf,nHei,nBlocks):
3204 def setup(self,dataOut,groupList,path,file,filec):
3205 self.__dataReady = False
3198 self.__dataReady = False
3206 # new
3199 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
3207 self.nChannels = dataOut.nChannels
3200 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
3208 self.channels = dataOut.channelList
3209 self.nHeights = dataOut.heightList.size
3210 self.heights = dataOut.heightList
3211 self.nProf = dataOut.nProfiles
3212 self.nIncohInt = dataOut.nIncohInt
3213 self.absc = dataOut.abscissaList[:-1]
3214
3215
3216 #To be inserted as a parameter
3217 try:
3218 self.groupArray = numpy.array(groupList)#groupArray = numpy.array([[0,1],[2,3]])
3219 except:
3220 print("Please insert groupList. Example (0,1),(2,3) format multilist")
3221 dataOut.groupList = self.groupArray
3222 self.crosspairs = dataOut.groupList
3223 self.nPairs = len(self.crosspairs)
3224 self.nGroups = self.groupArray.shape[0]
3225
3226 #List of possible combinations
3227
3228 self.listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2)
3229 self.indCross = numpy.zeros(len(list(self.listComb)), dtype = 'int')
3230
3231 #Parameters Array
3232 dataOut.data_param = None
3233 dataOut.data_paramC = None
3234 dataOut.clean_num_aver = None
3235 dataOut.coh_num_aver = None
3236 dataOut.tmp_spectra_i = None
3237 dataOut.tmp_cspectra_i = None
3238 dataOut.tmp_spectra_c = None
3239 dataOut.tmp_cspectra_c = None
3240 dataOut.index = None
3241
3242 if path != None:
3243 sys.path.append(path)
3244 self.library = importlib.import_module(file)
3245 if filec != None:
3246 self.weightf = importlib.import_module(filec)
3247
3248 #Set constants
3249 self.constants = self.library.setConstants(dataOut)
3250 dataOut.constants = self.constants
3251 self.M = dataOut.normFactor
3252 self.N = dataOut.nFFTPoints
3253 self.ippSeconds = dataOut.ippSeconds
3254 self.K = dataOut.nIncohInt
3255 self.pairsArray = numpy.array(dataOut.pairsList)
3256 self.snrth = 20
3257
3258 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):
3201 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):
3259
3202
3260 if (nicoh is None): nicoh = 1
3203 if (nicoh is None): nicoh = 1
@@ -3318,8 +3261,8 class SpectralFitting(Operation):
3318 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
3261 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
3319 snr = (spec2.mean()-n0)/n0
3262 snr = (spec2.mean()-n0)/n0
3320
3263
3321 if (snr < 1.e-20) :
3264 # if (snr < 1.e-20) :
3322 snr = 1.e-20
3265 # snr = 1.e-20
3323
3266
3324 vec_power[ind] = power
3267 vec_power[ind] = power
3325 vec_fd[ind] = fd
3268 vec_fd[ind] = fd
@@ -3329,8 +3272,10 class SpectralFitting(Operation):
3329 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
3272 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
3330 return moments
3273 return moments
3331
3274
3275 #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs):
3332 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
3276 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
3333
3277
3278 #import matplotlib.pyplot as plt
3334 nProf = dataOut.nProfiles
3279 nProf = dataOut.nProfiles
3335 heights = dataOut.heightList
3280 heights = dataOut.heightList
3336 nHei = len(heights)
3281 nHei = len(heights)
@@ -3363,10 +3308,12 class SpectralFitting(Operation):
3363 s_n1 = power[pair[1],:]/noise[pair[1]]
3308 s_n1 = power[pair[1],:]/noise[pair[1]]
3364 valid1 =(s_n0>=snr_th).nonzero()
3309 valid1 =(s_n0>=snr_th).nonzero()
3365 valid2 = (s_n1>=snr_th).nonzero()
3310 valid2 = (s_n1>=snr_th).nonzero()
3311
3366 valid1 = numpy.array(valid1[0])
3312 valid1 = numpy.array(valid1[0])
3367 valid2 = numpy.array(valid2[0])
3313 valid2 = numpy.array(valid2[0])
3368 valid = valid1
3314 valid = valid1
3369 for iv in range(len(valid2)):
3315 for iv in range(len(valid2)):
3316
3370 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3317 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3371 if len(indv[0]) == 0 :
3318 if len(indv[0]) == 0 :
3372 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
3319 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
@@ -3374,13 +3321,16 class SpectralFitting(Operation):
3374 my_coh_aver[pair[0],valid]=1
3321 my_coh_aver[pair[0],valid]=1
3375 my_coh_aver[pair[1],valid]=1
3322 my_coh_aver[pair[1],valid]=1
3376 # si la coherencia es mayor a la coherencia threshold los datos se toman
3323 # si la coherencia es mayor a la coherencia threshold los datos se toman
3324
3377 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)))
3325 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)))
3326
3378 for ih in range(len(hei_th)):
3327 for ih in range(len(hei_th)):
3379 hvalid = (heights>hei_th[ih]).nonzero()
3328 hvalid = (heights>hei_th[ih]).nonzero()
3380 hvalid = hvalid[0]
3329 hvalid = hvalid[0]
3381 if len(hvalid)>0:
3330 if len(hvalid)>0:
3382 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
3331 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
3383 valid = valid[0]
3332 valid = valid[0]
3333
3384 if len(valid)>0:
3334 if len(valid)>0:
3385 my_coh_aver[pair[0],hvalid[valid]] =1
3335 my_coh_aver[pair[0],hvalid[valid]] =1
3386 my_coh_aver[pair[1],hvalid[valid]] =1
3336 my_coh_aver[pair[1],hvalid[valid]] =1
@@ -3415,6 +3365,7 class SpectralFitting(Operation):
3415 valid1 = numpy.array(valid1[0])
3365 valid1 = numpy.array(valid1[0])
3416 valid2 = numpy.array(valid2[0])
3366 valid2 = numpy.array(valid2[0])
3417 incoh_echoes = valid1
3367 incoh_echoes = valid1
3368
3418 for iv in range(len(valid2)):
3369 for iv in range(len(valid2)):
3419
3370
3420 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3371 indv = numpy.array((valid1 == valid2[iv]).nonzero())
@@ -3433,11 +3384,12 class SpectralFitting(Operation):
3433 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
3384 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
3434 incoh_aver[pair[0],incoh_echoes]=1
3385 incoh_aver[pair[0],incoh_echoes]=1
3435 incoh_aver[pair[1],incoh_echoes]=1
3386 incoh_aver[pair[1],incoh_echoes]=1
3436 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
3437
3387
3388 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
3438
3389
3439 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3390 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3440
3391
3392 #import matplotlib.pyplot as plt
3441 nProf = dataOut.nProfiles
3393 nProf = dataOut.nProfiles
3442 heights = dataOut.heightList
3394 heights = dataOut.heightList
3443 nHei = len(heights)
3395 nHei = len(heights)
@@ -3448,6 +3400,7 class SpectralFitting(Operation):
3448
3400
3449 absc = dataOut.abscissaList[:-1]
3401 absc = dataOut.abscissaList[:-1]
3450 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
3402 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
3403
3451 clean_coh_spectra = spectra.copy()
3404 clean_coh_spectra = spectra.copy()
3452 clean_coh_cspectra = cspectra.copy()
3405 clean_coh_cspectra = cspectra.copy()
3453 clean_coh_aver = coh_aver.copy()
3406 clean_coh_aver = coh_aver.copy()
@@ -3462,7 +3415,9 class SpectralFitting(Operation):
3462 if clean_coh_echoes == 1 :
3415 if clean_coh_echoes == 1 :
3463 for ind in range(nChan):
3416 for ind in range(nChan):
3464 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
3417 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
3418
3465 spwd = data_param[:,3]
3419 spwd = data_param[:,3]
3420
3466 # 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
3421 # 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
3467 # para obtener spwd
3422 # para obtener spwd
3468 for ic in range(nPairs):
3423 for ic in range(nPairs):
@@ -3507,17 +3462,21 class SpectralFitting(Operation):
3507 nPairs = len(crosspairs)
3462 nPairs = len(crosspairs)
3508 hval=(heights >= min_hei).nonzero()
3463 hval=(heights >= min_hei).nonzero()
3509 ih=hval[0]
3464 ih=hval[0]
3465
3466
3510 for ih in range(hval[0][0],nHei):
3467 for ih in range(hval[0][0],nHei):
3511 for ifreq in range(nProf):
3468 for ifreq in range(nProf):
3512 for ii in range(n_funct):
3469 for ii in range(n_funct):
3513
3470
3514 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3471 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3472
3515 val = (numpy.isfinite(func2clean)==True).nonzero()
3473 val = (numpy.isfinite(func2clean)==True).nonzero()
3516 if len(val)>0:
3474 if len(val)>0:
3517 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3475 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3518 if min_val <= -40 : min_val = -40
3476 if min_val <= -40 : min_val = -40
3519 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
3477 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
3520 if max_val >= 200 : max_val = 200
3478 if max_val >= 200 : max_val = 200
3479
3521 step = 1
3480 step = 1
3522 #Getting bins and the histogram
3481 #Getting bins and the histogram
3523 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
3482 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
@@ -3549,8 +3508,29 class SpectralFitting(Operation):
3549 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
3508 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
3550 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
3509 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
3551 cspectra[noval,ii,ifreq,ih] = numpy.nan
3510 cspectra[noval,ii,ifreq,ih] = numpy.nan
3511 #
3512 #no sale es para savedrifts >2
3513 ''' channels = dataOut.channelList
3514 cross_pairs = dataOut.groupList
3552
3515
3516 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
3517 vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
3518 vcross = numpy.concatenate((vcross0,vcross1),axis=None)
3519
3520 #Getting coherent echoes which are removed.
3521 if len(novall) > 0:
3522 #val_spc[novall,ii,ifreq,ih] = 1
3523 val_spc[ii,ifreq,ih,novall] = 1
3524 if len(vcross) > 0:
3525 val_cspc[vcross,ifreq,ih,novall] = 1
3526
3527 #Removing coherent from ISR data.
3528 spectra[ii,ifreq,ih,noval] = numpy.nan
3529 if len(vcross) > 0:
3530 cspectra[vcross,ifreq,ih,noval] = numpy.nan
3531 '''
3553 #Getting average of the spectra and cross-spectra from incoherent echoes.
3532 #Getting average of the spectra and cross-spectra from incoherent echoes.
3533
3554 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
3534 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
3555 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3535 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3556 for ih in range(nHei):
3536 for ih in range(nHei):
@@ -3560,6 +3540,7 class SpectralFitting(Operation):
3560 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3540 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3561 if len(valid[0]) >0 :
3541 if len(valid[0]) >0 :
3562 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3542 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3543
3563 for icr in range(nPairs):
3544 for icr in range(nPairs):
3564 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3545 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3565 valid = (numpy.isfinite(tmp)==True).nonzero()
3546 valid = (numpy.isfinite(tmp)==True).nonzero()
@@ -3582,11 +3563,25 class SpectralFitting(Operation):
3582 for k in range(nHei):
3563 for k in range(nHei):
3583 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
3564 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
3584 val_cspc[:,i,j,k] = 0.0
3565 val_cspc[:,i,j,k] = 0.0
3566 # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHei*nChan))
3567 # if numpy.isfinite(val_spectra)==str(True):
3568 # noval = (val_spectra<1).nonzero()
3569 # if len(noval) > 0:
3570 # val_spc[:,noval] = 0.0
3571 # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHei))
3572
3573 #val_cspc = numpy.reshape(val_spc, (149,nChan*nHei*nProf))
3574 #if numpy.isfinite(val_cspectra)==str(True):
3575 # noval = (val_cspectra<1).nonzero()
3576 # if len(noval) > 0:
3577 # val_cspc[:,noval] = 0.0
3578 # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHei))
3585
3579
3586 tmp_sat_spectra = spectra.copy()
3580 tmp_sat_spectra = spectra.copy()
3587 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
3581 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
3588 tmp_sat_cspectra = cspectra.copy()
3582 tmp_sat_cspectra = cspectra.copy()
3589 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
3583 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
3584
3590 val = (val_spc > 0).nonzero()
3585 val = (val_spc > 0).nonzero()
3591 if len(val[0]) > 0:
3586 if len(val[0]) > 0:
3592 tmp_sat_spectra[val] = in_sat_spectra[val]
3587 tmp_sat_spectra[val] = in_sat_spectra[val]
@@ -3611,8 +3606,8 class SpectralFitting(Operation):
3611 valid = (numpy.isfinite(tmp)).nonzero()
3606 valid = (numpy.isfinite(tmp)).nonzero()
3612 if len(valid[0]) > 0:
3607 if len(valid[0]) > 0:
3613 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3608 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3614 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3615
3609
3610 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3616 def REM_ISOLATED_POINTS(self,array,rth):
3611 def REM_ISOLATED_POINTS(self,array,rth):
3617 if rth == None : rth = 4
3612 if rth == None : rth = 4
3618 num_prof = len(array[0,:,0])
3613 num_prof = len(array[0,:,0])
@@ -3627,24 +3622,28 class SpectralFitting(Operation):
3627 indxs1 = (indxs1[0])
3622 indxs1 = (indxs1[0])
3628 indxs2 = indxs2[0]
3623 indxs2 = indxs2[0]
3629 indxs = None
3624 indxs = None
3625
3630 for iv in range(len(indxs2)):
3626 for iv in range(len(indxs2)):
3631 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3627 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3632 if len(indv[0]) > 0 :
3628 if len(indv[0]) > 0 :
3633 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3629 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3630
3634 indxs = indxs[1:]
3631 indxs = indxs[1:]
3635 if len(indxs) < 4 :
3632 if len(indxs) < 4 :
3636 array[ii,:,:] = 0.
3633 array[ii,:,:] = 0.
3637 return
3634 return
3638
3635
3639 xpos = numpy.mod(indxs ,num_hei)
3636 xpos = numpy.mod(indxs ,num_prof)
3640 ypos = (indxs / num_hei)
3637 ypos = (indxs / num_prof)
3641 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3638 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3642 xpos = xpos[sx]
3639 xpos = xpos[sx]
3643 ypos = ypos[sx]
3640 ypos = ypos[sx]
3641
3644 # *********************************** Cleaning isolated points **********************************
3642 # *********************************** Cleaning isolated points **********************************
3645 ic = 0
3643 ic = 0
3646 while True :
3644 while True :
3647 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3645 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3646
3648 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3647 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3649 no_coh2 = (r <= rth).nonzero()
3648 no_coh2 = (r <= rth).nonzero()
3650 no_coh1 = numpy.array(no_coh1[0])
3649 no_coh1 = numpy.array(no_coh1[0])
@@ -3676,8 +3675,8 class SpectralFitting(Operation):
3676
3675
3677 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3676 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3678 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3677 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3679 return array
3680
3678
3679 return array
3681 def moments(self,doppler,yarray,npoints):
3680 def moments(self,doppler,yarray,npoints):
3682 ytemp = yarray
3681 ytemp = yarray
3683 val = (ytemp > 0).nonzero()
3682 val = (ytemp > 0).nonzero()
@@ -3695,19 +3694,23 class SpectralFitting(Operation):
3695 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3694 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3696 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3695 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3697 return [fmom,numpy.sqrt(smom)]
3696 return [fmom,numpy.sqrt(smom)]
3698
3697 # **********************************************************************************************
3698 index = 0
3699 fint = 0
3700 buffer = 0
3701 buffer2 = 0
3702 buffer3 = 0
3699 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None):
3703 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None,taver=None,proc=None,nhei=None,nprofs=None,ipp=None,channelList=None):
3700 if not self.isConfig:
3701 self.setup(dataOut = dataOut,groupList=groupList,path=path,file=file,filec=filec)
3702 self.isConfig = True
3703
3704 if not numpy.any(proc):
3704 if not numpy.any(proc):
3705 if numpy.any(taver):
3705
3706 taver = int(taver)
3706 nChannels = dataOut.nChannels
3707 else :
3707 nHeights= dataOut.heightList.size
3708 taver = 5
3708 nProf = dataOut.nProfiles
3709 if numpy.any(taver): taver=int(taver)
3710 else : taver = 5
3709 tini=time.localtime(dataOut.utctime)
3711 tini=time.localtime(dataOut.utctime)
3710 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3712 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3713
3711 self.index = 0
3714 self.index = 0
3712 jspc = self.buffer
3715 jspc = self.buffer
3713 jcspc = self.buffer2
3716 jcspc = self.buffer2
@@ -3717,18 +3720,15 class SpectralFitting(Operation):
3717 self.buffer3 = dataOut.noise
3720 self.buffer3 = dataOut.noise
3718 self.fint = 1
3721 self.fint = 1
3719 if numpy.any(jspc) :
3722 if numpy.any(jspc) :
3720 jspc = numpy.reshape(jspc ,(int(len(jspc) / self.nChannels) , self.nChannels ,self.nProf,self.nHeights ))
3723 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
3721 jcspc = numpy.reshape(jcspc ,(int(len(jcspc) /int(self.nChannels/2)),int(self.nChannels/2),self.nProf,self.nHeights ))
3724 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
3722 jnoise = numpy.reshape(jnoise,(int(len(jnoise)/ self.nChannels) , self.nChannels))
3725 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
3723 else:
3726 else:
3724 dataOut.flagNoData = True
3727 dataOut.flagNoData = True
3725 return dataOut
3728 return dataOut
3726 else :
3729 else :
3727 if (tini.tm_min % taver) == 0 :
3730 if (tini.tm_min % taver) == 0 : self.fint = 1
3728 self.fint = 1
3731 else : self.fint = 0
3729 else :
3730 self.fint = 0
3731
3732 self.index += 1
3732 self.index += 1
3733 if numpy.any(self.buffer):
3733 if numpy.any(self.buffer):
3734 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
3734 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
@@ -3740,29 +3740,106 class SpectralFitting(Operation):
3740 self.buffer3 = dataOut.noise
3740 self.buffer3 = dataOut.noise
3741 dataOut.flagNoData = True
3741 dataOut.flagNoData = True
3742 return dataOut
3742 return dataOut
3743 if path != None:
3744 sys.path.append(path)
3745 self.library = importlib.import_module(file)
3746 if filec != None:
3747 self.weightf = importlib.import_module(filec)
3748 #self.weightf = importlib.import_module('weightfit')
3749
3750 #To be inserted as a parameter
3751 groupArray = numpy.array(groupList)
3752 #groupArray = numpy.array([[0,1],[2,3]])
3753 dataOut.groupList = groupArray
3754
3755 nGroups = groupArray.shape[0]
3756 nChannels = dataOut.nChannels
3757 nHeights = dataOut.heightList.size
3758
3759 #Parameters Array
3760 dataOut.data_param = None
3761 dataOut.data_paramC = None
3762 dataOut.clean_num_aver = None
3763 dataOut.coh_num_aver = None
3764 dataOut.tmp_spectra_i = None
3765 dataOut.tmp_cspectra_i = None
3766 dataOut.tmp_spectra_c = None
3767 dataOut.tmp_cspectra_c = None
3768 dataOut.sat_spectra = None
3769 dataOut.sat_cspectra = None
3770 dataOut.index = None
3743
3771
3744 jnoise = jnoise/self.N# creo que falta dividirlo entre N
3772 #Set constants
3773 constants = self.library.setConstants(dataOut)
3774 dataOut.constants = constants
3775 M = dataOut.normFactor
3776 N = dataOut.nFFTPoints
3777
3778 ippSeconds = dataOut.ippSeconds
3779 K = dataOut.nIncohInt
3780 pairsArray = numpy.array(dataOut.pairsList)
3781
3782 snrth= 15
3783 spectra = dataOut.data_spc
3784 cspectra = dataOut.data_cspc
3785 nProf = dataOut.nProfiles
3786 heights = dataOut.heightList
3787 nHei = len(heights)
3788 channels = dataOut.channelList
3789 nChan = len(channels)
3790 nIncohInt = dataOut.nIncohInt
3791 crosspairs = dataOut.groupList
3792 noise = dataOut.noise
3793 jnoise = jnoise/N
3745 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
3794 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
3795 #print("NOISE-> ", 10*numpy.log10(noise))
3796 power = numpy.sum(spectra, axis=1)
3797 nPairs = len(crosspairs)
3798 absc = dataOut.abscissaList[:-1]
3799 #print('para escribir h5 ',dataOut.paramInterval)
3800 if not self.isConfig:
3801 self.isConfig = True
3802
3746 index = tini.tm_hour*12+tini.tm_min/taver
3803 index = tini.tm_hour*12+tini.tm_min/taver
3747 dataOut.index= index
3804 dataOut.index= index
3748 jspc = jspc/self.N/self.N
3805 jspc = jspc/N/N
3749 jcspc = jcspc/self.N/self.N
3806 jcspc = jcspc/N/N
3750
3751 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3807 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3752 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3808 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3753 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3809 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
3754
3810 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,coh_th, hei_th)
3755 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, self.snrth,coh_th, hei_th)
3811 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3756 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(self.snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3757
3758 dataOut.data_spc = incoh_spectra
3812 dataOut.data_spc = incoh_spectra
3759 dataOut.data_cspc = incoh_cspectra
3813 dataOut.data_cspc = incoh_cspectra
3814 dataOut.sat_spectra = sat_spectra
3815 dataOut.sat_cspectra = sat_cspectra
3816 # dataOut.data_spc = tmp_spectra
3817 # dataOut.data_cspc = tmp_cspectra
3818
3760 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3819 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
3761 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3820 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
3821 # clean_num_aver = (numpy.zeros([nChan, nHei])+1)*len(jspc[:,0,0,0])
3822 # coh_num_aver = numpy.zeros([nChan, nHei])*0*len(jspc[:,0,0,0])
3762 dataOut.clean_num_aver = clean_num_aver
3823 dataOut.clean_num_aver = clean_num_aver
3763 dataOut.coh_num_aver = coh_num_aver
3824 dataOut.coh_num_aver = coh_num_aver
3825 dataOut.tmp_spectra_i = incoh_spectra
3826 dataOut.tmp_cspectra_i = incoh_cspectra
3827 dataOut.tmp_spectra_c = clean_coh_spectra
3828 dataOut.tmp_cspectra_c = clean_coh_cspectra
3829 #List of possible combinations
3830 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3831 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3764
3832
3833 if getSNR:
3834 listChannels = groupArray.reshape((groupArray.size))
3835 listChannels.sort()
3836 norm = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter #* jspc.shape[0]
3837 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels], norm=norm)
3765 else:
3838 else:
3839 if numpy.any(taver): taver=int(taver)
3840 else : taver = 5
3841 tini=time.localtime(dataOut.utctime)
3842 index = tini.tm_hour*12+tini.tm_min/taver
3766 clean_num_aver = dataOut.clean_num_aver
3843 clean_num_aver = dataOut.clean_num_aver
3767 coh_num_aver = dataOut.coh_num_aver
3844 coh_num_aver = dataOut.coh_num_aver
3768 dataOut.data_spc = dataOut.tmp_spectra_i
3845 dataOut.data_spc = dataOut.tmp_spectra_i
@@ -3776,6 +3853,7 class SpectralFitting(Operation):
3776 dataOut.data_param = None
3853 dataOut.data_param = None
3777 dataOut.data_paramC = None
3854 dataOut.data_paramC = None
3778 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
3855 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
3856 #dataOut.paramInterval = 2.0
3779 #M=600
3857 #M=600
3780 #N=200
3858 #N=200
3781 dataOut.flagDecodeData=True
3859 dataOut.flagDecodeData=True
@@ -3785,66 +3863,83 class SpectralFitting(Operation):
3785 dataOut.nIncohInt= int(dataOut.nIncohInt)
3863 dataOut.nIncohInt= int(dataOut.nIncohInt)
3786 dataOut.nProfiles = int(dataOut.nProfiles)
3864 dataOut.nProfiles = int(dataOut.nProfiles)
3787 dataOut.nCohInt = int(dataOut.nCohInt)
3865 dataOut.nCohInt = int(dataOut.nCohInt)
3866 #print('sale',dataOut.nProfiles,dataOut.nHeights)
3788 #dataOut.nFFTPoints=nprofs
3867 #dataOut.nFFTPoints=nprofs
3789 #dataOut.normFactor = nprofs
3868 #dataOut.normFactor = nprofs
3790 dataOut.channelList = channelList
3869 dataOut.channelList = channelList
3870 nChan = len(channelList)
3791 #dataOut.ippFactor=1
3871 #dataOut.ippFactor=1
3792 #ipp = ipp/150*1.e-3
3872 #ipp = ipp/150*1.e-3
3793 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
3873 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
3794 #dataOut.ippSeconds=ipp
3874 #dataOut.ippSeconds=ipp
3795 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
3875 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
3876 #print('sale 2',dataOut.ippSeconds,M,N)
3877 # print('Empieza procesamiento offline')
3796 if path != None:
3878 if path != None:
3797 sys.path.append(path)
3879 sys.path.append(path)
3798 self.library = importlib.import_module(file)
3880 self.library = importlib.import_module(file)
3799 constants = self.library.setConstants(dataOut)
3881 constants = self.library.setConstants(dataOut)
3800 constants['M'] = M
3882 constants['M'] = M
3801 dataOut.constants = constants
3883 dataOut.constants = constants
3884 if filec != None:
3885 self.weightf = importlib.import_module(filec)
3802
3886
3887 groupArray = numpy.array(groupList)
3888 dataOut.groupList = groupArray
3889 nGroups = groupArray.shape[0]
3803 #List of possible combinations
3890 #List of possible combinations
3804 listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2)
3891 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3805 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3892 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3806 if dataOut.data_paramC is None:
3893 if dataOut.data_paramC is None:
3807 dataOut.data_paramC = numpy.zeros((self.nGroups*4, self.nHeights ,2))*numpy.nan
3894 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
3808 for i in range(self.nGroups):
3895 dataOut.data_snr1_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
3809 coord = self.groupArray[i,:]
3896 # dataOut.smooth_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
3897 for i in range(nGroups):
3898 coord = groupArray[i,:]
3810 #Input data array
3899 #Input data array
3811 data = dataOut.data_spc[coord,:,:]/(self.M*self.N)
3900 data = dataOut.data_spc[coord,:,:]/(M*N)
3812 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3901 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3813
3902
3814 #Cross Spectra data array for Covariance Matrixes
3903 #Cross Spectra data array for Covariance Matrixes
3815 ind = 0
3904 ind = 0
3816 for pairs in listComb:
3905 for pairs in listComb:
3817 pairsSel = numpy.array([coord[x],coord[y]])
3906 pairsSel = numpy.array([coord[x],coord[y]])
3818 indCross[ind] = int(numpy.where(numpy.all(self.pairsArray == pairsSel, axis = 1))[0][0])
3907 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
3819 ind += 1
3908 ind += 1
3820 dataCross = dataOut.data_cspc[indCross,:,:]/(self.M*self.N)
3909 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
3821 dataCross = dataCross**2
3910 dataCross = dataCross**2
3822 nhei = self.nHeights
3911 nhei = nHeights
3823 poweri = numpy.sum(dataOut.data_spc[:,1:self.nProf-0,:],axis=1)/clean_num_aver[:,:]
3912 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
3824
3825 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3913 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
3826 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(self.nProf-1)
3914 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
3827 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(self.nProf-1)
3915 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
3828 n0 = n0i
3916 n0 = n0i
3829 n1= n1i
3917 n1= n1i
3830 my_noises[2*i+0] = n0
3918 my_noises[2*i+0] = n0
3831 my_noises[2*i+1] = n1
3919 my_noises[2*i+1] = n1
3832 snrth = -15.0 # -4 -16 -25
3920 snrth = -13.0 # -4 -16 -25
3833 snrth = 10**(snrth/10.0)
3921 snrth = 10**(snrth/10.0)
3834 jvelr = numpy.zeros(self.nHeights, dtype = 'float')
3922 jvelr = numpy.zeros(nHeights, dtype = 'float')
3923 #snr0 = numpy.zeros(nHeights, dtype = 'float')
3924 #snr1 = numpy.zeros(nHeights, dtype = 'float')
3835 hvalid = [0]
3925 hvalid = [0]
3836 coh2 = abs(dataOut.data_cspc[i,1:self.nProf,:])**2/(dataOut.data_spc[0+i*2,1:self.nProf-0,:]*dataOut.data_spc[1+i*2,1:self.nProf-0,:])
3837
3926
3838 for h in range(self.nHeights):
3927 coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
3928
3929 for h in range(nHeights):
3839 smooth = clean_num_aver[i+1,h]
3930 smooth = clean_num_aver[i+1,h]
3840 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3931 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3841 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3932 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3842 signal0 = signalpn0-n0
3933 signal0 = signalpn0-n0
3843 signal1 = signalpn1-n1
3934 signal1 = signalpn1-n1
3844 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3935 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3845 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3936 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3937 #jmax0 = MAX(signal0,maxp0)
3938 #jmax1 = MAX(signal1,maxp1)
3846 gamma = coh2[:,h]
3939 gamma = coh2[:,h]
3940
3847 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3941 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3942
3848 if len(indxs) >0:
3943 if len(indxs) >0:
3849 if numpy.nanmean(gamma) > 0.07:
3944 if numpy.nanmean(gamma) > 0.07:
3850 maxp0 = numpy.argmax(signal0*gamma)
3945 maxp0 = numpy.argmax(signal0*gamma)
@@ -3853,22 +3948,24 class SpectralFitting(Operation):
3853 else:
3948 else:
3854 maxp0 = numpy.argmax(signal0)
3949 maxp0 = numpy.argmax(signal0)
3855 maxp1 = numpy.argmax(signal1)
3950 maxp1 = numpy.argmax(signal1)
3856 jvelr[h] = (self.absc[maxp0]+self.absc[maxp1])/2.
3951 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
3857 else: jvelr[h] = self.absc[0]
3952 else: jvelr[h] = absc[0]
3858 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3953 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3859 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3954 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3860
3955
3861 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3956 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3862 else: fd0 = numpy.nan
3957 else: fd0 = numpy.nan
3863 for h in range(self.nHeights):
3958 #print(fd0)
3959 for h in range(nHeights):
3864 d = data[:,h]
3960 d = data[:,h]
3865 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
3961 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
3866 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3962 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3867 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3963 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3868 signal0 = signalpn0-n0
3964 signal0 = signalpn0-n0
3869 signal1 = signalpn1-n1
3965 signal1 = signalpn1-n1
3870 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3966 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3871 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3967 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3968
3872 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3969 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3873 #Covariance Matrix
3970 #Covariance Matrix
3874 D = numpy.diag(d**2)
3971 D = numpy.diag(d**2)
@@ -3886,6 +3983,8 class SpectralFitting(Operation):
3886 ind += 1
3983 ind += 1
3887 diagD = numpy.zeros(256)
3984 diagD = numpy.zeros(256)
3888
3985
3986 #Dinv=numpy.linalg.inv(D)
3987 #L=numpy.linalg.cholesky(Dinv)
3889 try:
3988 try:
3890 Dinv=numpy.linalg.inv(D)
3989 Dinv=numpy.linalg.inv(D)
3891 L=numpy.linalg.cholesky(Dinv)
3990 L=numpy.linalg.cholesky(Dinv)
@@ -3895,58 +3994,84 class SpectralFitting(Operation):
3895 LT=L.T
3994 LT=L.T
3896
3995
3897 dp = numpy.dot(LT,d)
3996 dp = numpy.dot(LT,d)
3997
3898 #Initial values
3998 #Initial values
3899 data_spc = dataOut.data_spc[coord,:,h]
3999 data_spc = dataOut.data_spc[coord,:,h]
3900 w = data_spc/data_spc
4000 w = data_spc/data_spc
3901 if filec != None:
4001 if filec != None:
3902 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
4002 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
4003
3903 if (h>6) and (error1[3]<25):
4004 if (h>6) and (error1[3]<25):
3904 p0 = dataOut.data_param[i,:,h-1].copy()
4005 p0 = dataOut.data_param[i,:,h-1].copy()
4006 #print('usa anterior')
3905 else:
4007 else:
3906 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, self.constants))# sin el i(data_spc, constants, i)
4008 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
3907 p0[3] = fd0
4009 p0[3] = fd0
3908
3909 if filec != None:
4010 if filec != None:
3910 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
4011 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
3911
4012
4013 #if index == 175 and i==1 and h>=27 and h<=35: p0[3]=30
4014 #if h >= 6 and i==1 and h<= 10: print(p0)
3912 try:
4015 try:
3913 #Least Squares
4016 #Least Squares
3914 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,self.constants),full_output=True)
4017 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
3915 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
4018 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3916 #Chi square error
4019 #Chi square error
3917 error0 = numpy.sum(infodict['fvec']**2)/(2*self.N)
4020 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
3918 #Error with Jacobian
4021 #Error with Jacobian
3919 error1 = self.library.errorFunction(minp,self.constants,LT)
4022 error1 = self.library.errorFunction(minp,constants,LT)
4023 #if h >= 0 and h<= 10 and i ==0: print(p0,minp,error1)
4024 #if i>=0 and h>=0: print(index,h,minp[3])
4025 # print self.__residFunction(p0,dp,LT, constants)
4026 # print infodict['fvec']
4027 # print self.__residFunction(minp,dp,LT,constants)
3920
4028
3921 except:
4029 except:
3922 minp = p0*numpy.nan
4030 minp = p0*numpy.nan
3923 error0 = numpy.nan
4031 error0 = numpy.nan
3924 error1 = p0*numpy.nan
4032 error1 = p0*numpy.nan
4033 # s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0))
4034 # covp = covp*s_sq
4035 # error = []
4036 # for ip in range(len(minp)):
4037 # try:
4038 # error.append(numpy.absolute(covp[ip][ip])**0.5)
4039 # except:
4040 # error.append( 0.00 )
4041 #if i==1 and h==11 and index == 139: print(p0, minp,data_spc)
3925 else :
4042 else :
3926 data_spc = dataOut.data_spc[coord,:,h]
4043 data_spc = dataOut.data_spc[coord,:,h]
3927 p0 = numpy.array(self.library.initialValuesFunction(data_spc, self.constants))
4044 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
3928 minp = p0*numpy.nan
4045 minp = p0*numpy.nan
3929 error0 = numpy.nan
4046 error0 = numpy.nan
3930 error1 = p0*numpy.nan
4047 error1 = p0*numpy.nan
4048
3931 if dataOut.data_param is None:
4049 if dataOut.data_param is None:
3932 dataOut.data_param = numpy.zeros((self.nGroups, p0.size, self.nHeights ))*numpy.nan
4050 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
3933 dataOut.data_error = numpy.zeros((self.nGroups, p0.size + 1, self.nHeights ))*numpy.nan
4051 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
4052
3934 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
4053 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3935 dataOut.data_param[i,:,h] = minp
4054 dataOut.data_param[i,:,h] = minp
3936
4055 dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0
3937 for ht in range(self.nHeights-1) :
4056 dataOut.data_snr1_i[i*2+1,h] = numpy.sum(signalpn1/(nProf-1))/n1
4057 #dataOut.smooth_i[i*2,h] = clean_num_aver[i+1,h]
4058 #print(fd0,dataOut.data_param[i,3,h])
4059 #print(fd0,dataOut.data_param[i,3,:])
4060 for ht in range(nHeights-1) :
3938 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
4061 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3939 dataOut.data_paramC[4*i,ht,1] = smooth
4062 dataOut.data_paramC[4*i,ht,1] = smooth
3940 signalpn0 = (clean_coh_spectra[i*2 ,1:(self.nProf-0),ht])/smooth #coh_spectra
4063 signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
3941 signalpn1 = (clean_coh_spectra[i*2+1,1:(self.nProf-0),ht])/smooth
4064 signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
4065
3942 val0 = (signalpn0 > 0).nonzero()
4066 val0 = (signalpn0 > 0).nonzero()
3943 val0 = val0[0]
4067 val0 = val0[0]
3944 if len(val0) == 0 : val0_npoints = self.nProf
4068
4069 if len(val0) == 0 : val0_npoints = nProf
3945 else : val0_npoints = len(val0)
4070 else : val0_npoints = len(val0)
3946
4071
3947 val1 = (signalpn1 > 0).nonzero()
4072 val1 = (signalpn1 > 0).nonzero()
3948 val1 = val1[0]
4073 val1 = val1[0]
3949 if len(val1) == 0 : val1_npoints = self.nProf
4074 if len(val1) == 0 : val1_npoints = nProf
3950 else : val1_npoints = len(val1)
4075 else : val1_npoints = len(val1)
3951
4076
3952 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
4077 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
@@ -3960,31 +4085,62 class SpectralFitting(Operation):
3960 vali = (signal1 < 0).nonzero()
4085 vali = (signal1 < 0).nonzero()
3961 vali = vali[0]
4086 vali = vali[0]
3962 if len(vali) > 0 : signal1[vali] = 0
4087 if len(vali) > 0 : signal1[vali] = 0
3963 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
4088 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3964 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
4089 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3965 doppler = self.absc[1:]
4090 doppler = absc[1:]
3966 if snr0 >= snrth and snr1 >= snrth and smooth :
4091 if snr0 >= snrth and snr1 >= snrth and smooth :
3967 signalpn0_n0 = signalpn0
4092 signalpn0_n0 = signalpn0
3968 signalpn0_n0[val0] = signalpn0[val0] - n0
4093 signalpn0_n0[val0] = signalpn0[val0] - n0
3969 mom0 = self.moments(doppler,signalpn0-n0,self.nProf)
4094 mom0 = self.moments(doppler,signalpn0-n0,nProf)
4095
3970 signalpn1_n1 = signalpn1
4096 signalpn1_n1 = signalpn1
3971 signalpn1_n1[val1] = signalpn1[val1] - n1
4097 signalpn1_n1[val1] = signalpn1[val1] - n1
3972 mom1 = self.moments(doppler,signalpn1_n1,self.nProf)
4098 mom1 = self.moments(doppler,signalpn1_n1,nProf)
3973 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
4099 dataOut.data_paramC[2+4*i,ht,0] = (mom0[0]+mom1[0])/2.
3974 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
4100 dataOut.data_paramC[3+4*i,ht,0] = (mom0[1]+mom1[1])/2.
4101 #dataOut.data_snr1_c[i*2,ht] = numpy.sum(signalpn0/(nProf-1))/n0
4102 #dataOut.data_snr1_c[i*2+1,ht] = numpy.sum(signalpn1/(nProf-1))/n1
3975 dataOut.data_spc = jspectra
4103 dataOut.data_spc = jspectra
3976 dataOut.spc_noise = my_noises*self.nProf*self.M
4104 dataOut.spc_noise = my_noises*nProf*M
3977 if numpy.any(proc): dataOut.spc_noise = my_noises*self.nProf*self.M
4105
3978 if getSNR:
4106 if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
3979 listChannels = self.groupArray.reshape((self.groupArray.size))
4107 if 0:
4108 listChannels = groupArray.reshape((groupArray.size))
3980 listChannels.sort()
4109 listChannels.sort()
3981 # TEST
4110 norm = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
3982 noise_C = numpy.zeros(self.nChannels)
4111 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels], norm=norm)
3983 noise_C = dataOut.getNoise()
4112 #print(dataOut.data_snr1_i)
3984 #print("noise_C",noise_C)
4113 # Adding coherent echoes from possible satellites.
3985 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:],noise_C/(600.0*1.15))# PRUEBA *nProf*M
4114 #sat_spectra = numpy.zeros((nChan,nProf,nHei), dtype=float)
3986 #dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise_C[listChannels])# PRUEBA *nProf*M
4115 #sat_spectra = sat_spectra[*,*,anal_header.channels]
3987 dataOut.flagNoData = False
4116 isat_spectra = numpy.zeros([2,int(nChan/2),nProf,nhei], dtype=float)
4117
4118 sat_fits = numpy.zeros([4,nhei], dtype=float)
4119 noises = my_noises/nProf
4120 #nchan2 = int(nChan/2)
4121 for beam in range(int(nChan/2)-0) :
4122 n0 = noises[2*beam]
4123 n1 = noises[2*beam+1]
4124 isat_spectra[0:2,beam,:,:] = dataOut.sat_spectra[2*beam +0:2*beam+2 ,:,:]
4125
4126 for ht in range(nhei-1) :
4127 signalpn0 = isat_spectra[0,beam,:,ht]
4128 signalpn0 = numpy.reshape(signalpn0,nProf)
4129 signalpn1 = isat_spectra[1,beam,:,ht]
4130 signalpn1 = numpy.reshape(signalpn1,nProf)
4131
4132 cval0 = len((signalpn0 > 0).nonzero()[0])
4133 if cval0 == 0 : val0_npoints = nProf
4134 else: val0_npoints = cval0
4135
4136 cval1 = len((signalpn1 > 0).nonzero()[0])
4137 if cval1 == 0 : val1_npoints = nProf
4138 else: val1_npoints = cval1
4139
4140 sat_fits[0+2*beam,ht] = numpy.sum(signalpn0/(val0_npoints*nProf))/n0
4141 sat_fits[1+2*beam,ht] = numpy.sum(signalpn1/(val1_npoints*nProf))/n1
4142
4143 dataOut.sat_fits = sat_fits
3988 return dataOut
4144 return dataOut
3989
4145
3990 def __residFunction(self, p, dp, LT, constants):
4146 def __residFunction(self, p, dp, LT, constants):
@@ -3993,10 +4149,16 class SpectralFitting(Operation):
3993 fmp=numpy.dot(LT,fm)
4149 fmp=numpy.dot(LT,fm)
3994 return dp-fmp
4150 return dp-fmp
3995
4151
3996 def __getSNR(self, z, noise):
4152 def __getSNR(self, z, noise, norm=1):
3997
4153
4154 # normFactor = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * pwcode * dataOut.windowOfFilter
3998 avg = numpy.average(z, axis=1)
4155 avg = numpy.average(z, axis=1)
4156 #noise /= norm
3999 SNR = (avg.T-noise)/noise
4157 SNR = (avg.T-noise)/noise
4158 # SNR = avg.T/noise
4159 # print("Noise: ", 10*numpy.log10(noise))
4160 # print("SNR: ", SNR)
4161 # print("SNR: ", 10*numpy.log10(SNR))
4000 SNR = SNR.T
4162 SNR = SNR.T
4001 return SNR
4163 return SNR
4002
4164
@@ -4007,7 +4169,6 class SpectralFitting(Operation):
4007 fmp=numpy.dot(LT,fm)
4169 fmp=numpy.dot(LT,fm)
4008 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
4170 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
4009 return chisq
4171 return chisq
4010
4011 class WindProfiler(Operation):
4172 class WindProfiler(Operation):
4012
4173
4013 __isConfig = False
4174 __isConfig = False
@@ -4671,6 +4832,7 class EWDriftsEstimation(Operation):
4671 y1=y1[vali]
4832 y1=y1[vali]
4672 x = x[vali]
4833 x = x[vali]
4673 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
4834 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
4835
4674 x1 = heiRang1
4836 x1 = heiRang1
4675 y11 = f1(x1)
4837 y11 = f1(x1)
4676 y2 = SNR[i,:]
4838 y2 = SNR[i,:]
@@ -4678,6 +4840,7 class EWDriftsEstimation(Operation):
4678 vali= (y2 != -1).nonzero()
4840 vali= (y2 != -1).nonzero()
4679 y2 = y2[vali]
4841 y2 = y2[vali]
4680 x = x[vali]
4842 x = x[vali]
4843
4681 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
4844 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
4682 y21 = f2(x1)
4845 y21 = f2(x1)
4683
4846
@@ -4686,17 +4849,38 class EWDriftsEstimation(Operation):
4686
4849
4687 return heiRang1, velRadial1, SNR1
4850 return heiRang1, velRadial1, SNR1
4688
4851
4689 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4690
4852
4853
4854 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4691 dataOut.lat = -11.95
4855 dataOut.lat = -11.95
4692 dataOut.lon = -76.87
4856 dataOut.lon = -76.87
4857 dataOut.spcst = 0.00666
4858 dataOut.pl = 0.0003
4859 dataOut.cbadn = 3
4860 dataOut.inttms = 300
4861 dataOut.azw = -115.687
4862 dataOut.elw = 86.1095
4863 dataOut.aze = 130.052
4864 dataOut.ele = 87.6558
4865 dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor)
4866 dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor)
4867 dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor)
4868 dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor)
4869
4693 heiRang = dataOut.heightList
4870 heiRang = dataOut.heightList
4694 velRadial = dataOut.data_param[:,3,:]
4871 velRadial = dataOut.data_param[:,3,:]
4695 velRadialm = dataOut.data_param[:,2:4,:]*-1
4872 velRadialm = dataOut.data_param[:,2:4,:]*-1
4873
4696 rbufc=dataOut.data_paramC[:,:,0]
4874 rbufc=dataOut.data_paramC[:,:,0]
4697 ebufc=dataOut.data_paramC[:,:,1]
4875 ebufc=dataOut.data_paramC[:,:,1]
4698 SNR = dataOut.data_snr
4876 SNR = dataOut.data_snr1_i
4877 rbufi = dataOut.data_snr1_i
4699 velRerr = dataOut.data_error[:,4,:]
4878 velRerr = dataOut.data_error[:,4,:]
4879 range1 = dataOut.heightList
4880 nhei = len(range1)
4881
4882 sat_fits = dataOut.sat_fits
4883
4700 channels = dataOut.channelList
4884 channels = dataOut.channelList
4701 nChan = len(channels)
4885 nChan = len(channels)
4702 my_nbeams = nChan/2
4886 my_nbeams = nChan/2
@@ -4705,18 +4889,48 class EWDriftsEstimation(Operation):
4705 else :
4889 else :
4706 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
4890 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
4707 dataOut.moments=moments
4891 dataOut.moments=moments
4892 #Incoherent
4893 smooth_w = dataOut.clean_num_aver[0,:]
4894 chisq_w = dataOut.data_error[0,0,:]
4895 p_w0 = rbufi[0,:]
4896 p_w1 = rbufi[1,:]
4708 # Coherent
4897 # Coherent
4709 smooth_wC = ebufc[0,:]
4898 smooth_wC = ebufc[0,:]
4710 p_w0C = rbufc[0,:]
4899 p_w0C = rbufc[0,:]
4711 p_w1C = rbufc[1,:]
4900 p_w1C = rbufc[1,:]
4712 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
4901 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
4713 t_wC = rbufc[3,:]
4902 t_wC = rbufc[3,:]
4903 val = (numpy.isfinite(p_w0)==False).nonzero()
4904 p_w0[val]=0
4905 val = (numpy.isfinite(p_w1)==False).nonzero()
4906 p_w1[val]=0
4907 val = (numpy.isfinite(p_w0C)==False).nonzero()
4908 p_w0C[val]=0
4909 val = (numpy.isfinite(p_w1C)==False).nonzero()
4910 p_w1C[val]=0
4911 val = (numpy.isfinite(smooth_w)==False).nonzero()
4912 smooth_w[val]=0
4913 val = (numpy.isfinite(smooth_wC)==False).nonzero()
4914 smooth_wC[val]=0
4915
4916 #p_w0 = (p_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC)
4917 #p_w1 = (p_w1*smooth_w+p_w1C*smooth_wC)/(smooth_w+smooth_wC)
4918
4919 if len(sat_fits) >0 :
4920 p_w0C = p_w0C + sat_fits[0,:]
4921 p_w1C = p_w1C + sat_fits[1,:]
4922
4714 if my_nbeams == 1:
4923 if my_nbeams == 1:
4715 w = velRadial[0,:]
4924 w = velRadial[0,:]
4716 winds = velRadial.copy()
4925 winds = velRadial.copy()
4717 w_err = velRerr[0,:]
4926 w_err = velRerr[0,:]
4718 snr1 = 10*numpy.log10(SNR[0])
4927 u = w*numpy.nan
4928 u_err = w_err*numpy.nan
4929 p_e0 = p_w0*numpy.nan
4930 p_e1 = p_w1*numpy.nan
4931 #snr1 = 10*numpy.log10(SNR[0])
4719 if my_nbeams == 2:
4932 if my_nbeams == 2:
4933
4720 zenith = numpy.array(zenith)
4934 zenith = numpy.array(zenith)
4721 zenith -= zenithCorrection
4935 zenith -= zenithCorrection
4722 zenith *= numpy.pi/180
4936 zenith *= numpy.pi/180
@@ -4734,24 +4948,72 class EWDriftsEstimation(Operation):
4734 w_e = velRadial1[1,:]
4948 w_e = velRadial1[1,:]
4735 w_w_err = velRerr[0,:]
4949 w_w_err = velRerr[0,:]
4736 w_e_err = velRerr[1,:]
4950 w_e_err = velRerr[1,:]
4951 smooth_e = dataOut.clean_num_aver[2,:]
4952 chisq_e = dataOut.data_error[1,0,:]
4953 p_e0 = rbufi[2,:]
4954 p_e1 = rbufi[3,:]
4737
4955
4738 val = (numpy.isfinite(w_w)==False).nonzero()
4956 tini=time.localtime(dataOut.utctime)
4739 val = val[0]
4957 #print(tini[3],tini[4])
4740 bad = val
4958 #val = (numpy.isfinite(w_w)==False).nonzero()
4741 if len(bad) > 0 :
4959 #val = val[0]
4742 w_w[bad] = w_wC[bad]
4960 #bad = val
4743 w_w_err[bad]= numpy.nan
4961 #if len(bad) > 0 :
4962 # w_w[bad] = w_wC[bad]
4963 # w_w_err[bad]= numpy.nan
4964 if tini[3] >= 6 and tini[3] < 18 :
4965 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
4966 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
4967 else:
4968 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
4969 w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp)
4970 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
4971 w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp)
4972 w_w = w_wtmp
4973 w_w_err = w_w_errtmp
4974
4975 #if my_nbeams == 2:
4744 smooth_eC=ebufc[4,:]
4976 smooth_eC=ebufc[4,:]
4745 p_e0C = rbufc[4,:]
4977 p_e0C = rbufc[4,:]
4746 p_e1C = rbufc[5,:]
4978 p_e1C = rbufc[5,:]
4747 w_eC = rbufc[6,:]*-1
4979 w_eC = rbufc[6,:]*-1
4748 t_eC = rbufc[7,:]
4980 t_eC = rbufc[7,:]
4749 val = (numpy.isfinite(w_e)==False).nonzero()
4981
4750 val = val[0]
4982 val = (numpy.isfinite(p_e0)==False).nonzero()
4751 bad = val
4983 p_e0[val]=0
4752 if len(bad) > 0 :
4984 val = (numpy.isfinite(p_e1)==False).nonzero()
4753 w_e[bad] = w_eC[bad]
4985 p_e1[val]=0
4754 w_e_err[bad]= numpy.nan
4986 val = (numpy.isfinite(p_e0C)==False).nonzero()
4987 p_e0C[val]=0
4988 val = (numpy.isfinite(p_e1C)==False).nonzero()
4989 p_e1C[val]=0
4990 val = (numpy.isfinite(smooth_e)==False).nonzero()
4991 smooth_e[val]=0
4992 val = (numpy.isfinite(smooth_eC)==False).nonzero()
4993 smooth_eC[val]=0
4994 #p_e0 = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC)
4995 #p_e1 = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC)
4996
4997 if len(sat_fits) >0 :
4998 p_e0C = p_e0C + sat_fits[2,:]
4999 p_e1C = p_e1C + sat_fits[3,:]
5000
5001 #val = (numpy.isfinite(w_e)==False).nonzero()
5002 #val = val[0]
5003 #bad = val
5004 #if len(bad) > 0 :
5005 # w_e[bad] = w_eC[bad]
5006 # w_e_err[bad]= numpy.nan
5007 if tini[3] >= 6 and tini[3] < 18 :
5008 w_etmp = numpy.where(numpy.isfinite(w_eC)==True,w_eC,w_e)
5009 w_e_errtmp = numpy.where(numpy.isfinite(w_eC)==True,numpy.nan,w_e_err)
5010 else:
5011 w_etmp = numpy.where(numpy.isfinite(w_eC)==True,w_eC,w_e)
5012 w_etmp = numpy.where(range1 > 200,w_e,w_etmp)
5013 w_e_errtmp = numpy.where(numpy.isfinite(w_eC)==True,numpy.nan,w_e_err)
5014 w_e_errtmp = numpy.where(range1 > 200,w_e_err,w_e_errtmp)
5015 w_e = w_etmp
5016 w_e_err = w_e_errtmp
4755
5017
4756 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
5018 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
4757 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
5019 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
@@ -4760,13 +5022,22 class EWDriftsEstimation(Operation):
4760 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))
5022 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))
4761
5023
4762 winds = numpy.vstack((w,u))
5024 winds = numpy.vstack((w,u))
4763
5025 #winds = numpy.vstack((w,u,w_err,u_err))
4764 dataOut.heightList = heiRang1
5026 dataOut.heightList = heiRang1
4765 snr1 = 10*numpy.log10(SNR1[0])
5027 #snr1 = 10*numpy.log10(SNR1[0])
4766 dataOut.data_output = winds
5028 dataOut.data_output = winds
4767 #snr1 = 10*numpy.log10(SNR1[0])# estaba comentado
5029 range1 = dataOut.heightList
4768 dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
5030 nhei = len(range1)
4769 #print("data_snr1",dataOut.data_snr1)
5031 #print('alt ',range1*numpy.sin(86.1*numpy.pi/180))
5032 #print(numpy.min([dataOut.eldir7,dataOut.eldir8]))
5033 galt = range1*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.)
5034 dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e))
5035 #snr1 = 10*numpy.log10(SNR1[0])
5036 #print(min(snr1), max(snr1))
5037 snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1))
5038 snr1db = 10*numpy.log10(snr1[0])
5039 #dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
5040 dataOut.data_snr1 = numpy.reshape(snr1db,(1,snr1db.shape[0]))
4770 dataOut.utctimeInit = dataOut.utctime
5041 dataOut.utctimeInit = dataOut.utctime
4771 dataOut.outputInterval = dataOut.timeInterval
5042 dataOut.outputInterval = dataOut.timeInterval
4772
5043
@@ -4774,8 +5045,6 class EWDriftsEstimation(Operation):
4774 jrange = 450 #900 para HA drifts
5045 jrange = 450 #900 para HA drifts
4775 deltah = 15.0 #dataOut.spacing(0) 25 HAD
5046 deltah = 15.0 #dataOut.spacing(0) 25 HAD
4776 h0 = 0.0 #dataOut.first_height(0)
5047 h0 = 0.0 #dataOut.first_height(0)
4777 heights = dataOut.heightList
4778 nhei = len(heights)
4779
5048
4780 range1 = numpy.arange(nhei) * deltah + h0
5049 range1 = numpy.arange(nhei) * deltah + h0
4781 jhei = (range1 >= hei_aver0).nonzero()
5050 jhei = (range1 >= hei_aver0).nonzero()
@@ -4840,19 +5109,96 class EWDriftsEstimation(Operation):
4840 uA = uA[6*h_avgs:nhei_avg-0]
5109 uA = uA[6*h_avgs:nhei_avg-0]
4841 uA_err = uA_err[6*h_avgs:nhei_avg-0]
5110 uA_err = uA_err[6*h_avgs:nhei_avg-0]
4842 dataOut.drifts_avg = numpy.vstack((wA,uA))
5111 dataOut.drifts_avg = numpy.vstack((wA,uA))
5112
4843 if my_nbeams == 1: dataOut.drifts_avg = wA
5113 if my_nbeams == 1: dataOut.drifts_avg = wA
5114 #deltahavg= wA*0.0+deltah
5115 dataOut.range = range1
5116 galtavg = range_aver*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.)
5117 dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h))
5118
5119 #print('comparando dim de avg ',wA.shape,deltahavg.shape,range_aver.shape)
4844 tini=time.localtime(dataOut.utctime)
5120 tini=time.localtime(dataOut.utctime)
4845 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
5121 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
4846 nfile = fileDrifts+'/jro'+datefile+'drifts.txt'
5122 nfile = fileDrifts+'/jro'+datefile+'drifts_sch3.txt'
5123 #nfile = '/home/pcondor/Database/ewdriftsschain2019/jro'+datefile+'drifts_sch3.txt'
5124 #print(len(dataOut.drifts_avg),dataOut.drifts_avg.shape)
4847 f1 = open(nfile,'a')
5125 f1 = open(nfile,'a')
4848 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
5126 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
4849 driftavgstr=str(dataOut.drifts_avg)
5127 driftavgstr=str(dataOut.drifts_avg)
4850 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
5128 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
4851 numpy.savetxt(f1,dataOut.drifts_avg,fmt='%10.2f')
5129 numpy.savetxt(f1,numpy.reshape(range_aver,(1,len(range_aver))) ,fmt='%10.2f')
5130 #numpy.savetxt(f1,numpy.reshape(dataOut.drifts_avg,(7,len(dataOut.drifts_avg))) ,fmt='%10.2f')
5131 numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f')
4852 f1.close()
5132 f1.close()
4853
5133
5134 swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt'
5135 f1 = open(swfile,'a')
5136 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
5137 numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f')
5138 numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f')
5139 f1.close()
5140 dataOut.heightListtmp = dataOut.heightList
5141 '''
5142 one = {'range':'range','gdlatr': 'lat', 'gdlonr': 'lon', 'inttms': 'paramInterval'} #reader gdlatr-->lat only 1D
5143
5144 two = {
5145 'gdalt': 'heightList', #<----- nmonics
5146 'VIPN': ('params', 0),
5147 'dvipn': ('params', 1),
5148 'vipe': ('params', 2),
5149 'dvipe': ('params', 3),
5150 'PACWL': ('params', 4),
5151 'pbcwl': ('params', 5),
5152 'pccel': ('params', 6),
5153 'pdcel': ('params', 7)
5154 } #writer
5155
5156 #f=open('/home/roberto/moder_test.txt','r')
5157 #file_contents=f.read()
5158
5159 ind = ['gdalt']
5160
5161 meta = {
5162 'kinst': 10, #instrument code
5163 'kindat': 1910, #type of data
5164 'catalog': {
5165 'principleInvestigator': 'Danny ScipiΓ³n',
5166 'expPurpose': 'Drifts'#,
5167 # 'sciRemarks': file_contents
5168 },
5169 'header': {
5170 'analyst': 'D. Hysell'
5171 }
5172 }
5173 print('empieza h5 madrigal')
5174 try:
5175 h5mad=MADWriter(dataOut, fileDrifts, one, ind, two, meta, format='hdf5')
5176 except:
5177 print("Error in MADWriter")
5178 print(h5mad)
5179 '''
5180 return dataOut
5181 class setHeightDrifts(Operation):
5182
5183 def __init__(self):
5184 Operation.__init__(self)
5185 def run(self, dataOut):
5186 #print('h inicial ',dataOut.heightList,dataOut.heightListtmp)
5187 dataOut.heightList = dataOut.heightListtmp
5188 #print('regresa H ',dataOut.heightList)
5189 return dataOut
5190 class setHeightDriftsavg(Operation):
5191
5192 def __init__(self):
5193 Operation.__init__(self)
5194 def run(self, dataOut):
5195 #print('h inicial ',dataOut.heightList)
5196 dataOut.heightList = dataOut.params_avg[4]
5197 #print('cambia H ',dataOut.params_avg[4],dataOut.heightList)
4854 return dataOut
5198 return dataOut
4855
5199
5200
5201
4856 #--------------- Non Specular Meteor ----------------
5202 #--------------- Non Specular Meteor ----------------
4857
5203
4858 class NonSpecularMeteorDetection(Operation):
5204 class NonSpecularMeteorDetection(Operation):
General Comments 0
You need to be logged in to leave comments. Login now