##// 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 459 powerdB = numpy.squeeze(powerdB)
460 460
461 461 return powerdB
462 @property
463 def data_pow(self):
464 return self.getPower()
462 465
463 466 @property
464 467 def timeInterval(self):
Show file before | Show file after | Hide comments
@@ -206,7 +206,7 class GenericRTIPlot(Plot):
206 206 data = {
207 207 'param' : numpy.concatenate([getattr(dataOut, attr) for attr in self.attr_data], axis=0)
208 208 }
209
209
210 210 meta = {}
211 211
212 212 return data, meta
Show file before | Show file after | Hide comments
@@ -65,22 +65,26 class SpectraPlot(Plot):
65 65 self.update_list(dataOut)
66 66 data = {}
67 67 meta = {}
68
69 68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
70 noise = 10*numpy.log10(dataOut.getNoise()/norm)
71 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
72 for ch in range(dataOut.nChannels):
73 if hasattr(dataOut.normFactor,'ndim'):
74 if dataOut.normFactor.ndim > 1:
75 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
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:
73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74
75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
76 for ch in range(dataOut.nChannels):
77 if hasattr(dataOut.normFactor,'ndim'):
78 if dataOut.normFactor.ndim > 1:
79 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
76 80
81 else:
82 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
77 83 else:
78 84 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
79 else:
80 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
81 85
82 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
83 spc = 10*numpy.log10(z)
86 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
87 spc = 10*numpy.log10(z)
84 88
85 89 data['spc'] = spc
86 90 data['rti'] = spc.mean(axis=1)
@@ -725,7 +729,7 class RTIPlot(Plot):
725 729 )
726 730 if self.showprofile:
727 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 733 if "noise" in self.data:
730 734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
731 735 color="k", linestyle="dashed", lw=1)[0]
@@ -737,7 +741,7 class RTIPlot(Plot):
737 741 cmap=plt.get_cmap(self.colormap)
738 742 )
739 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 745 if "noise" in self.data:
742 746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
743 747 color="k", linestyle="dashed", lw=1)[0]
@@ -893,12 +897,15 class NoisePlot(Plot):
893 897 self.titles = ['Noise']
894 898 self.colorbar = False
895 899 self.plots_adjust.update({'right': 0.85 })
900 self.titles = ['Noise Plot']
896 901
897 902 def update(self, dataOut):
898 903
899 904 data = {}
900 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 909 meta['yrange'] = numpy.array([])
903 910
904 911 return data, meta
Show file before | Show file after | Hide comments
@@ -445,6 +445,7 class HDFWriter(Operation):
445 445 path = None
446 446 setFile = None
447 447 fp = None
448 ds = None
448 449 firsttime = True
449 450 #Configurations
450 451 blocksPerFile = None
@@ -490,6 +491,8 class HDFWriter(Operation):
490 491 if self.metadataList is None:
491 492 self.metadataList = self.dataOut.metadata_list
492 493
494 self.metadataList = list(set(self.metadataList))
495
493 496 tableList = []
494 497 dsList = []
495 498
@@ -504,7 +507,7 class HDFWriter(Operation):
504 507
505 508 if dataAux is None:
506 509 continue
507 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
510 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float_)):
508 511 dsDict['nDim'] = 0
509 512 else:
510 513 dsDict['nDim'] = len(dataAux.shape)
@@ -514,6 +517,7 class HDFWriter(Operation):
514 517
515 518 dsList.append(dsDict)
516 519
520 self.blockIndex = 0
517 521 self.dsList = dsList
518 522 self.currentDay = self.dataOut.datatime.date()
519 523
@@ -763,7 +767,7 class HDFWriter(Operation):
763 767 self.ds = dtsets
764 768 self.data = data
765 769 self.firsttime = True
766 self.blockIndex = 0
770
767 771 return
768 772
769 773 def putData(self):
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (840 lines changed) Show them Hide them
@@ -24,7 +24,7 import warnings
24 24 from numpy import NaN
25 25 from scipy.optimize.optimize import OptimizeWarning
26 26 warnings.filterwarnings('ignore')
27
27 import json
28 28 import os
29 29 import csv
30 30 from scipy import signal
@@ -216,7 +216,7 class ParametersProc(ProcessingUnit):
216 216 self.__updateObjFromInput()
217 217 self.dataOut.utctimeInit = self.dataIn.utctime
218 218 self.dataOut.paramInterval = self.dataIn.timeInterval
219
219
220 220 return
221 221
222 222
@@ -3176,7 +3176,6 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
3176 3176 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
3177 3177 return y
3178 3178
3179
3180 3179 class SpectralFitting(Operation):
3181 3180 '''
3182 3181 Function GetMoments()
@@ -3189,72 +3188,16 class SpectralFitting(Operation):
3189 3188 __dataReady = False
3190 3189 bloques = None
3191 3190 bloque0 = None
3192 index = 0
3193 fint = 0
3194 buffer = 0
3195 buffer2 = 0
3196 buffer3 = 0
3197 3191
3198 3192 def __init__(self):
3199 3193 Operation.__init__(self)
3200 3194 self.i=0
3201 3195 self.isConfig = False
3202 3196
3203
3204 def setup(self,dataOut,groupList,path,file,filec):
3197 def setup(self,nChan,nProf,nHei,nBlocks):
3205 3198 self.__dataReady = False
3206 # new
3207 self.nChannels = dataOut.nChannels
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
3199 self.bloques = numpy.zeros([2, nProf, nHei,nBlocks], dtype= complex)
3200 self.bloque0 = numpy.zeros([nChan, nProf, nHei, nBlocks])
3258 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 3203 if (nicoh is None): nicoh = 1
@@ -3318,19 +3261,21 class SpectralFitting(Operation):
3318 3261 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
3319 3262 snr = (spec2.mean()-n0)/n0
3320 3263
3321 if (snr < 1.e-20) :
3322 snr = 1.e-20
3264 # if (snr < 1.e-20) :
3265 # snr = 1.e-20
3323 3266
3324 3267 vec_power[ind] = power
3325 3268 vec_fd[ind] = fd
3326 3269 vec_w[ind] = w
3327 3270 vec_snr[ind] = snr
3328
3271
3329 3272 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
3330 3273 return moments
3331 3274
3275 #def __DiffCoherent(self,snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise, crosspairs):
3332 3276 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
3333
3277
3278 #import matplotlib.pyplot as plt
3334 3279 nProf = dataOut.nProfiles
3335 3280 heights = dataOut.heightList
3336 3281 nHei = len(heights)
@@ -3363,10 +3308,12 class SpectralFitting(Operation):
3363 3308 s_n1 = power[pair[1],:]/noise[pair[1]]
3364 3309 valid1 =(s_n0>=snr_th).nonzero()
3365 3310 valid2 = (s_n1>=snr_th).nonzero()
3311
3366 3312 valid1 = numpy.array(valid1[0])
3367 3313 valid2 = numpy.array(valid2[0])
3368 3314 valid = valid1
3369 3315 for iv in range(len(valid2)):
3316
3370 3317 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3371 3318 if len(indv[0]) == 0 :
3372 3319 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
@@ -3374,17 +3321,20 class SpectralFitting(Operation):
3374 3321 my_coh_aver[pair[0],valid]=1
3375 3322 my_coh_aver[pair[1],valid]=1
3376 3323 # si la coherencia es mayor a la coherencia threshold los datos se toman
3324
3377 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 3327 for ih in range(len(hei_th)):
3379 3328 hvalid = (heights>hei_th[ih]).nonzero()
3380 3329 hvalid = hvalid[0]
3381 3330 if len(hvalid)>0:
3382 3331 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
3383 3332 valid = valid[0]
3333
3384 3334 if len(valid)>0:
3385 3335 my_coh_aver[pair[0],hvalid[valid]] =1
3386 3336 my_coh_aver[pair[1],hvalid[valid]] =1
3387
3337
3388 3338 coh_echoes = (my_coh_aver[pair[0],:] == 1).nonzero()
3389 3339 incoh_echoes = (my_coh_aver[pair[0],:] != 1).nonzero()
3390 3340 incoh_echoes = incoh_echoes[0]
@@ -3404,9 +3354,9 class SpectralFitting(Operation):
3404 3354 valid1 = numpy.array(valid1[0])
3405 3355 valid2 = numpy.array(valid2[0])
3406 3356 valid = valid1
3407
3357
3408 3358 for iv in range(len(valid2)):
3409
3359
3410 3360 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3411 3361 if len(indv[0]) == 0 :
3412 3362 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
@@ -3415,8 +3365,9 class SpectralFitting(Operation):
3415 3365 valid1 = numpy.array(valid1[0])
3416 3366 valid2 = numpy.array(valid2[0])
3417 3367 incoh_echoes = valid1
3368
3418 3369 for iv in range(len(valid2)):
3419
3370
3420 3371 indv = numpy.array((valid1 == valid2[iv]).nonzero())
3421 3372 if len(indv[0]) == 0 :
3422 3373 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
@@ -3433,11 +3384,12 class SpectralFitting(Operation):
3433 3384 incoh_cspectra[ic,:,incoh_echoes] = cspectra[ic,:,incoh_echoes]
3434 3385 incoh_aver[pair[0],incoh_echoes]=1
3435 3386 incoh_aver[pair[1],incoh_echoes]=1
3387
3436 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
3437
3438
3389
3439 3390 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
3440 3391
3392 #import matplotlib.pyplot as plt
3441 3393 nProf = dataOut.nProfiles
3442 3394 heights = dataOut.heightList
3443 3395 nHei = len(heights)
@@ -3448,6 +3400,7 class SpectralFitting(Operation):
3448 3400
3449 3401 absc = dataOut.abscissaList[:-1]
3450 3402 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
3403
3451 3404 clean_coh_spectra = spectra.copy()
3452 3405 clean_coh_cspectra = cspectra.copy()
3453 3406 clean_coh_aver = coh_aver.copy()
@@ -3462,7 +3415,9 class SpectralFitting(Operation):
3462 3415 if clean_coh_echoes == 1 :
3463 3416 for ind in range(nChan):
3464 3417 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
3418
3465 3419 spwd = data_param[:,3]
3420
3466 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 3422 # para obtener spwd
3468 3423 for ic in range(nPairs):
@@ -3487,9 +3442,9 class SpectralFitting(Operation):
3487 3442 clean_coh_aver[pair,ih] = 0
3488 3443
3489 3444 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
3490
3445
3491 3446 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
3492
3447
3493 3448 rfunc = cspectra.copy()
3494 3449 n_funct = len(rfunc[0,:,0,0])
3495 3450 val_spc = spectra*0.0
@@ -3507,19 +3462,23 class SpectralFitting(Operation):
3507 3462 nPairs = len(crosspairs)
3508 3463 hval=(heights >= min_hei).nonzero()
3509 3464 ih=hval[0]
3465
3466
3510 3467 for ih in range(hval[0][0],nHei):
3511 3468 for ifreq in range(nProf):
3512 3469 for ii in range(n_funct):
3513 3470
3514 3471 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
3472
3515 3473 val = (numpy.isfinite(func2clean)==True).nonzero()
3516 3474 if len(val)>0:
3517 3475 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
3518 3476 if min_val <= -40 : min_val = -40
3519 3477 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
3520 3478 if max_val >= 200 : max_val = 200
3479
3521 3480 step = 1
3522 #Getting bins and the histogram
3481 #Getting bins and the histogram
3523 3482 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
3524 3483 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
3525 3484 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
@@ -3532,7 +3491,7 class SpectralFitting(Operation):
3532 3491 except:
3533 3492 mode = mean
3534 3493 stdv = sigma
3535
3494
3536 3495 #Removing echoes greater than mode + 3*stdv
3537 3496 factor_stdv = 2.5
3538 3497 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
@@ -3541,31 +3500,53 class SpectralFitting(Operation):
3541 3500 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
3542 3501 cross_pairs = crosspairs[ii]
3543 3502 #Getting coherent echoes which are removed.
3544 if len(novall[0]) > 0:
3503 if len(novall[0]) > 0:
3545 3504 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
3546 3505 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
3547 val_cspc[novall[0],ii,ifreq,ih] = 1
3506 val_cspc[novall[0],ii,ifreq,ih] = 1
3548 3507 #Removing coherent from ISR data
3549 3508 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
3550 3509 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
3551 3510 cspectra[noval,ii,ifreq,ih] = numpy.nan
3552
3511 #
3512 #no sale es para savedrifts >2
3513 ''' channels = dataOut.channelList
3514 cross_pairs = dataOut.groupList
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 3532 #Getting average of the spectra and cross-spectra from incoherent echoes.
3533
3554 3534 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
3555 3535 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
3556 3536 for ih in range(nHei):
3557 3537 for ifreq in range(nProf):
3558 3538 for ich in range(nChan):
3559 3539 tmp = spectra[:,ich,ifreq,ih]
3560 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3540 valid = (numpy.isfinite(tmp[:])==True).nonzero()
3561 3541 if len(valid[0]) >0 :
3562 3542 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3543
3563 3544 for icr in range(nPairs):
3564 3545 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
3565 3546 valid = (numpy.isfinite(tmp)==True).nonzero()
3566 3547 if len(valid[0]) > 0:
3567 3548 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3568 #Removing fake coherent echoes (at least 4 points around the point)
3549 #Removing fake coherent echoes (at least 4 points around the point)
3569 3550 val_spectra = numpy.sum(val_spc,0)
3570 3551 val_cspectra = numpy.sum(val_cspc,0)
3571 3552
@@ -3582,11 +3563,25 class SpectralFitting(Operation):
3582 3563 for k in range(nHei):
3583 3564 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
3584 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 3580 tmp_sat_spectra = spectra.copy()
3587 3581 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
3588 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 3585 val = (val_spc > 0).nonzero()
3591 3586 if len(val[0]) > 0:
3592 3587 tmp_sat_spectra[val] = in_sat_spectra[val]
@@ -3611,8 +3606,8 class SpectralFitting(Operation):
3611 3606 valid = (numpy.isfinite(tmp)).nonzero()
3612 3607 if len(valid[0]) > 0:
3613 3608 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
3609
3614 3610 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
3615
3616 3611 def REM_ISOLATED_POINTS(self,array,rth):
3617 3612 if rth == None : rth = 4
3618 3613 num_prof = len(array[0,:,0])
@@ -3620,31 +3615,35 class SpectralFitting(Operation):
3620 3615 n2d = len(array[:,0,0])
3621 3616
3622 3617 for ii in range(n2d) :
3623 tmp = array[ii,:,:]
3618 tmp = array[ii,:,:]
3624 3619 tmp = numpy.reshape(tmp,num_prof*num_hei)
3625 3620 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
3626 indxs2 = (tmp > 0).nonzero()
3621 indxs2 = (tmp > 0).nonzero()
3627 3622 indxs1 = (indxs1[0])
3628 indxs2 = indxs2[0]
3623 indxs2 = indxs2[0]
3629 3624 indxs = None
3625
3630 3626 for iv in range(len(indxs2)):
3631 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3627 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
3632 3628 if len(indv[0]) > 0 :
3633 3629 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
3630
3634 3631 indxs = indxs[1:]
3635 3632 if len(indxs) < 4 :
3636 3633 array[ii,:,:] = 0.
3637 3634 return
3638 3635
3639 xpos = numpy.mod(indxs ,num_hei)
3640 ypos = (indxs / num_hei)
3636 xpos = numpy.mod(indxs ,num_prof)
3637 ypos = (indxs / num_prof)
3641 3638 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
3642 3639 xpos = xpos[sx]
3643 3640 ypos = ypos[sx]
3644 # *********************************** Cleaning isolated points **********************************
3641
3642 # *********************************** Cleaning isolated points **********************************
3645 3643 ic = 0
3646 3644 while True :
3647 3645 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
3646
3648 3647 no_coh1 = (numpy.isfinite(r)==True).nonzero()
3649 3648 no_coh2 = (r <= rth).nonzero()
3650 3649 no_coh1 = numpy.array(no_coh1[0])
@@ -3676,10 +3675,10 class SpectralFitting(Operation):
3676 3675
3677 3676 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3678 3677 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3678
3679 3679 return array
3680
3681 3680 def moments(self,doppler,yarray,npoints):
3682 ytemp = yarray
3681 ytemp = yarray
3683 3682 val = (ytemp > 0).nonzero()
3684 3683 val = val[0]
3685 3684 if len(val) == 0 : val = range(npoints-1)
@@ -3695,42 +3694,43 class SpectralFitting(Operation):
3695 3694 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3696 3695 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3697 3696 return [fmom,numpy.sqrt(smom)]
3698
3697 # **********************************************************************************************
3698 index = 0
3699 fint = 0
3700 buffer = 0
3701 buffer2 = 0
3702 buffer3 = 0
3699 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):
3705 if numpy.any(taver):
3706 taver = int(taver)
3707 else :
3708 taver = 5
3709 tini = time.localtime(dataOut.utctime)
3710 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3711 self.index = 0
3712 jspc = self.buffer
3713 jcspc = self.buffer2
3714 jnoise = self.buffer3
3715 self.buffer = dataOut.data_spc
3704 if not numpy.any(proc):
3705
3706 nChannels = dataOut.nChannels
3707 nHeights= dataOut.heightList.size
3708 nProf = dataOut.nProfiles
3709 if numpy.any(taver): taver=int(taver)
3710 else : taver = 5
3711 tini=time.localtime(dataOut.utctime)
3712 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
3713
3714 self.index = 0
3715 jspc = self.buffer
3716 jcspc = self.buffer2
3717 jnoise = self.buffer3
3718 self.buffer = dataOut.data_spc
3716 3719 self.buffer2 = dataOut.data_cspc
3717 3720 self.buffer3 = dataOut.noise
3718 self.fint = 1
3721 self.fint = 1
3719 3722 if numpy.any(jspc) :
3720 jspc = numpy.reshape(jspc ,(int(len(jspc) / self.nChannels) , self.nChannels ,self.nProf,self.nHeights ))
3721 jcspc = numpy.reshape(jcspc ,(int(len(jcspc) /int(self.nChannels/2)),int(self.nChannels/2),self.nProf,self.nHeights ))
3722 jnoise = numpy.reshape(jnoise,(int(len(jnoise)/ self.nChannels) , self.nChannels))
3723 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
3724 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
3725 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
3723 3726 else:
3724 3727 dataOut.flagNoData = True
3725 3728 return dataOut
3726 3729 else :
3727 if (tini.tm_min % taver) == 0 :
3728 self.fint = 1
3729 else :
3730 self.fint = 0
3731
3730 if (tini.tm_min % taver) == 0 : self.fint = 1
3731 else : self.fint = 0
3732 3732 self.index += 1
3733 if numpy.any(self.buffer):
3733 if numpy.any(self.buffer):
3734 3734 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
3735 3735 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
3736 3736 self.buffer3 = numpy.concatenate((self.buffer3,dataOut.noise), axis=0)
@@ -3740,29 +3740,106 class SpectralFitting(Operation):
3740 3740 self.buffer3 = dataOut.noise
3741 3741 dataOut.flagNoData = True
3742 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
3743 3758
3744 jnoise = jnoise/self.N# creo que falta dividirlo entre N
3745 noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
3746 index = tini.tm_hour*12+tini.tm_min/taver
3747 dataOut.index = index
3748 jspc = jspc/self.N/self.N
3749 jcspc = jcspc/self.N/self.N
3750
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
3771
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
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
3803 index = tini.tm_hour*12+tini.tm_min/taver
3804 dataOut.index= index
3805 jspc = jspc/N/N
3806 jcspc = jcspc/N/N
3751 3807 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
3752 jspectra = tmp_spectra * len(jspc[:,0,0,0])
3753 jcspectra = tmp_cspectra * len(jspc[:,0,0,0])
3754
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)
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
3759 dataOut.data_cspc = incoh_cspectra
3760 clean_num_aver = incoh_aver * len(jspc[:,0,0,0])
3761 coh_num_aver = clean_coh_aver* len(jspc[:,0,0,0])
3762 dataOut.clean_num_aver = clean_num_aver
3763 dataOut.coh_num_aver = coh_num_aver
3764
3808 jspectra = tmp_spectra*len(jspc[:,0,0,0])
3809 jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
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)
3811 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
3812 dataOut.data_spc = incoh_spectra
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
3819 clean_num_aver = incoh_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])
3823 dataOut.clean_num_aver = clean_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')
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 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 3843 clean_num_aver = dataOut.clean_num_aver
3767 3844 coh_num_aver = dataOut.coh_num_aver
3768 3845 dataOut.data_spc = dataOut.tmp_spectra_i
@@ -3776,6 +3853,7 class SpectralFitting(Operation):
3776 3853 dataOut.data_param = None
3777 3854 dataOut.data_paramC = None
3778 3855 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
3856 #dataOut.paramInterval = 2.0
3779 3857 #M=600
3780 3858 #N=200
3781 3859 dataOut.flagDecodeData=True
@@ -3785,66 +3863,83 class SpectralFitting(Operation):
3785 3863 dataOut.nIncohInt= int(dataOut.nIncohInt)
3786 3864 dataOut.nProfiles = int(dataOut.nProfiles)
3787 3865 dataOut.nCohInt = int(dataOut.nCohInt)
3866 #print('sale',dataOut.nProfiles,dataOut.nHeights)
3788 3867 #dataOut.nFFTPoints=nprofs
3789 3868 #dataOut.normFactor = nprofs
3790 3869 dataOut.channelList = channelList
3870 nChan = len(channelList)
3791 3871 #dataOut.ippFactor=1
3792 3872 #ipp = ipp/150*1.e-3
3793 3873 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
3794 3874 #dataOut.ippSeconds=ipp
3795 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 3878 if path != None:
3797 3879 sys.path.append(path)
3798 3880 self.library = importlib.import_module(file)
3799 3881 constants = self.library.setConstants(dataOut)
3800 3882 constants['M'] = M
3801 3883 dataOut.constants = constants
3802
3803 #List of possible combinations
3804 listComb = itertools.combinations(numpy.arange(self.groupArray.shape[1]),2)
3884 if filec != None:
3885 self.weightf = importlib.import_module(filec)
3886
3887 groupArray = numpy.array(groupList)
3888 dataOut.groupList = groupArray
3889 nGroups = groupArray.shape[0]
3890 #List of possible combinations
3891 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
3805 3892 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
3806 3893 if dataOut.data_paramC is None:
3807 dataOut.data_paramC = numpy.zeros((self.nGroups*4, self.nHeights ,2))*numpy.nan
3808 for i in range(self.nGroups):
3809 coord = self.groupArray[i,:]
3894 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
3895 dataOut.data_snr1_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
3896 # dataOut.smooth_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
3897 for i in range(nGroups):
3898 coord = groupArray[i,:]
3810 3899 #Input data array
3811 data = dataOut.data_spc[coord,:,:]/(self.M*self.N)
3812 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3900 data = dataOut.data_spc[coord,:,:]/(M*N)
3901 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
3813 3902
3814 3903 #Cross Spectra data array for Covariance Matrixes
3815 3904 ind = 0
3816 3905 for pairs in listComb:
3817 pairsSel = numpy.array([coord[x],coord[y]])
3818 indCross[ind] = int(numpy.where(numpy.all(self.pairsArray == pairsSel, axis = 1))[0][0])
3819 ind += 1
3820 dataCross = dataOut.data_cspc[indCross,:,:]/(self.M*self.N)
3906 pairsSel = numpy.array([coord[x],coord[y]])
3907 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
3908 ind += 1
3909 dataCross = dataOut.data_cspc[indCross,:,:]/(M*N)
3821 3910 dataCross = dataCross**2
3822 nhei = self.nHeights
3823 poweri = numpy.sum(dataOut.data_spc[:,1:self.nProf-0,:],axis=1)/clean_num_aver[:,:]
3824
3911 nhei = nHeights
3912 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
3825 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)
3827 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(self.nProf-1)
3914 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
3915 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
3828 3916 n0 = n0i
3829 3917 n1= n1i
3830 3918 my_noises[2*i+0] = n0
3831 3919 my_noises[2*i+1] = n1
3832 snrth = -15.0 # -4 -16 -25
3920 snrth = -13.0 # -4 -16 -25
3833 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 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
3838 for h in range(self.nHeights):
3926
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 3930 smooth = clean_num_aver[i+1,h]
3840 signalpn0 = (dataOut.data_spc[i*2,1:(self.nProf-0),h])/smooth
3841 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3931 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3932 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3842 3933 signal0 = signalpn0-n0
3843 3934 signal1 = signalpn1-n1
3844 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3845 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3935 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3936 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3937 #jmax0 = MAX(signal0,maxp0)
3938 #jmax1 = MAX(signal1,maxp1)
3846 3939 gamma = coh2[:,h]
3940
3847 3941 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
3942
3848 3943 if len(indxs) >0:
3849 3944 if numpy.nanmean(gamma) > 0.07:
3850 3945 maxp0 = numpy.argmax(signal0*gamma)
@@ -3853,24 +3948,26 class SpectralFitting(Operation):
3853 3948 else:
3854 3949 maxp0 = numpy.argmax(signal0)
3855 3950 maxp1 = numpy.argmax(signal1)
3856 jvelr[h] = (self.absc[maxp0]+self.absc[maxp1])/2.
3857 else: jvelr[h] = self.absc[0]
3951 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
3952 else: jvelr[h] = absc[0]
3858 3953 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
3859 3954 #print(maxp0,absc[maxp0],snr0,jvelr[h])
3860 3955
3861 3956 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
3862 3957 else: fd0 = numpy.nan
3863 for h in range(self.nHeights):
3958 #print(fd0)
3959 for h in range(nHeights):
3864 3960 d = data[:,h]
3865 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
3867 signalpn1 = (dataOut.data_spc[i*2+1,1:(self.nProf-0),h])/smooth
3962 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
3963 signalpn1 = (dataOut.data_spc[i*2+1,1:(nProf-0),h])/smooth
3868 3964 signal0 = signalpn0-n0
3869 3965 signal1 = signalpn1-n1
3870 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3871 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3966 snr0 = numpy.sum(signal0/n0)/(nProf-1)
3967 snr1 = numpy.sum(signal1/n1)/(nProf-1)
3968
3872 3969 if snr0 > snrth and snr1 > snrth and clean_num_aver[i+1,h] > 0 :
3873 #Covariance Matrix
3970 #Covariance Matrix
3874 3971 D = numpy.diag(d**2)
3875 3972 ind = 0
3876 3973 for pairs in listComb:
@@ -3885,7 +3982,9 class SpectralFitting(Operation):
3885 3982 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
3886 3983 ind += 1
3887 3984 diagD = numpy.zeros(256)
3888
3985
3986 #Dinv=numpy.linalg.inv(D)
3987 #L=numpy.linalg.cholesky(Dinv)
3889 3988 try:
3890 3989 Dinv=numpy.linalg.inv(D)
3891 3990 L=numpy.linalg.cholesky(Dinv)
@@ -3895,58 +3994,84 class SpectralFitting(Operation):
3895 3994 LT=L.T
3896 3995
3897 3996 dp = numpy.dot(LT,d)
3898 #Initial values
3997
3998 #Initial values
3899 3999 data_spc = dataOut.data_spc[coord,:,h]
3900 4000 w = data_spc/data_spc
3901 4001 if filec != None:
3902 4002 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
3903 if (h>6)and(error1[3]<25):
4003
4004 if (h>6) and (error1[3]<25):
3904 4005 p0 = dataOut.data_param[i,:,h-1].copy()
4006 #print('usa anterior')
3905 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 4009 p0[3] = fd0
3908
3909 4010 if filec != None:
3910 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 4015 try:
3913 #Least Squares
3914 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,self.constants),full_output=True)
3915 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
3916 #Chi square error
3917 error0 = numpy.sum(infodict['fvec']**2)/(2*self.N)
3918 #Error with Jacobian
3919 error1 = self.library.errorFunction(minp,self.constants,LT)
4016 #Least Squares
4017 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
4018 #minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
4019 #Chi square error
4020 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
4021 #Error with Jacobian
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 4029 except:
3922 4030 minp = p0*numpy.nan
3923 4031 error0 = numpy.nan
3924 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 4042 else :
3926 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 4045 minp = p0*numpy.nan
3929 4046 error0 = numpy.nan
3930 error1 = p0*numpy.nan
4047 error1 = p0*numpy.nan
4048
3931 4049 if dataOut.data_param is None:
3932 dataOut.data_param = numpy.zeros((self.nGroups, p0.size, self.nHeights ))*numpy.nan
3933 dataOut.data_error = numpy.zeros((self.nGroups, p0.size + 1, self.nHeights ))*numpy.nan
4050 dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
4051 dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
4052
3934 4053 dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
3935 4054 dataOut.data_param[i,:,h] = minp
3936
3937 for ht in range(self.nHeights-1) :
4055 dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0
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 4061 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
3939 4062 dataOut.data_paramC[4*i,ht,1] = smooth
3940 signalpn0 = (clean_coh_spectra[i*2 ,1:(self.nProf-0),ht])/smooth #coh_spectra
3941 signalpn1 = (clean_coh_spectra[i*2+1,1:(self.nProf-0),ht])/smooth
4063 signalpn0 = (clean_coh_spectra[i*2 ,1:(nProf-0),ht])/smooth #coh_spectra
4064 signalpn1 = (clean_coh_spectra[i*2+1,1:(nProf-0),ht])/smooth
4065
3942 4066 val0 = (signalpn0 > 0).nonzero()
3943 4067 val0 = val0[0]
3944 if len(val0) == 0 : val0_npoints = self.nProf
4068
4069 if len(val0) == 0 : val0_npoints = nProf
3945 4070 else : val0_npoints = len(val0)
3946
4071
3947 4072 val1 = (signalpn1 > 0).nonzero()
3948 4073 val1 = val1[0]
3949 if len(val1) == 0 : val1_npoints = self.nProf
4074 if len(val1) == 0 : val1_npoints = nProf
3950 4075 else : val1_npoints = len(val1)
3951 4076
3952 4077 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
@@ -3960,31 +4085,62 class SpectralFitting(Operation):
3960 4085 vali = (signal1 < 0).nonzero()
3961 4086 vali = vali[0]
3962 4087 if len(vali) > 0 : signal1[vali] = 0
3963 snr0 = numpy.sum(signal0/n0)/(self.nProf-1)
3964 snr1 = numpy.sum(signal1/n1)/(self.nProf-1)
3965 doppler = self.absc[1:]
4088 snr0 = numpy.sum(signal0/n0)/(nProf-1)
4089 snr1 = numpy.sum(signal1/n1)/(nProf-1)
4090 doppler = absc[1:]
3966 4091 if snr0 >= snrth and snr1 >= snrth and smooth :
3967 4092 signalpn0_n0 = signalpn0
3968 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 4096 signalpn1_n1 = signalpn1
3971 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 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 4103 dataOut.data_spc = jspectra
3976 dataOut.spc_noise = my_noises*self.nProf*self.M
3977 if numpy.any(proc): dataOut.spc_noise = my_noises*self.nProf*self.M
3978 if getSNR:
3979 listChannels = self.groupArray.reshape((self.groupArray.size))
4104 dataOut.spc_noise = my_noises*nProf*M
4105
4106 if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
4107 if 0:
4108 listChannels = groupArray.reshape((groupArray.size))
3980 4109 listChannels.sort()
3981 # TEST
3982 noise_C = numpy.zeros(self.nChannels)
3983 noise_C = dataOut.getNoise()
3984 #print("noise_C",noise_C)
3985 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:],noise_C/(600.0*1.15))# PRUEBA *nProf*M
3986 #dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise_C[listChannels])# PRUEBA *nProf*M
3987 dataOut.flagNoData = False
4110 norm = dataOut.nProfiles * dataOut.nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
4111 dataOut.data_snr = self.__getSNR(dataOut.data_spc[listChannels,:,:], my_noises[listChannels], norm=norm)
4112 #print(dataOut.data_snr1_i)
4113 # Adding coherent echoes from possible satellites.
4114 #sat_spectra = numpy.zeros((nChan,nProf,nHei), dtype=float)
4115 #sat_spectra = sat_spectra[*,*,anal_header.channels]
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 4144 return dataOut
3989 4145
3990 4146 def __residFunction(self, p, dp, LT, constants):
@@ -3993,10 +4149,16 class SpectralFitting(Operation):
3993 4149 fmp=numpy.dot(LT,fm)
3994 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 4155 avg = numpy.average(z, axis=1)
4156 #noise /= norm
3999 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 4162 SNR = SNR.T
4001 4163 return SNR
4002 4164
@@ -4006,8 +4168,7 class SpectralFitting(Operation):
4006 4168 dp=numpy.dot(LT,d)
4007 4169 fmp=numpy.dot(LT,fm)
4008 4170 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
4009 return chisq
4010
4171 return chisq
4011 4172 class WindProfiler(Operation):
4012 4173
4013 4174 __isConfig = False
@@ -4655,7 +4816,7 class EWDriftsEstimation(Operation):
4655 4816 maxid = listPhi.index(max(listPhi))
4656 4817 minid = listPhi.index(min(listPhi))
4657 4818
4658 rango = list(range(len(phi)))
4819 rango = list(range(len(phi)))
4659 4820 heiRang1 = heiRang*math.cos(phi[maxid])
4660 4821 heiRangAux = heiRang*math.cos(phi[minid])
4661 4822 indOut = (heiRang1 < heiRangAux[0]).nonzero()
@@ -4671,13 +4832,15 class EWDriftsEstimation(Operation):
4671 4832 y1=y1[vali]
4672 4833 x = x[vali]
4673 4834 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
4835
4674 4836 x1 = heiRang1
4675 4837 y11 = f1(x1)
4676 y2 = SNR[i,:]
4838 y2 = SNR[i,:]
4677 4839 x = heiRang*math.cos(phi[i])
4678 4840 vali= (y2 != -1).nonzero()
4679 4841 y2 = y2[vali]
4680 4842 x = x[vali]
4843
4681 4844 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
4682 4845 y21 = f2(x1)
4683 4846
@@ -4686,17 +4849,38 class EWDriftsEstimation(Operation):
4686 4849
4687 4850 return heiRang1, velRadial1, SNR1
4688 4851
4852
4853
4689 4854 def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
4690
4691 dataOut.lat=-11.95
4692 dataOut.lon=-76.87
4855 dataOut.lat = -11.95
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 4870 heiRang = dataOut.heightList
4694 4871 velRadial = dataOut.data_param[:,3,:]
4695 4872 velRadialm = dataOut.data_param[:,2:4,:]*-1
4873
4696 4874 rbufc=dataOut.data_paramC[:,:,0]
4697 4875 ebufc=dataOut.data_paramC[:,:,1]
4698 SNR = dataOut.data_snr
4876 SNR = dataOut.data_snr1_i
4877 rbufi = dataOut.data_snr1_i
4699 4878 velRerr = dataOut.data_error[:,4,:]
4879 range1 = dataOut.heightList
4880 nhei = len(range1)
4881
4882 sat_fits = dataOut.sat_fits
4883
4700 4884 channels = dataOut.channelList
4701 4885 nChan = len(channels)
4702 4886 my_nbeams = nChan/2
@@ -4705,18 +4889,48 class EWDriftsEstimation(Operation):
4705 4889 else :
4706 4890 moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
4707 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 4897 # Coherent
4709 4898 smooth_wC = ebufc[0,:]
4710 4899 p_w0C = rbufc[0,:]
4711 4900 p_w1C = rbufc[1,:]
4712 4901 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
4713 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 4923 if my_nbeams == 1:
4715 4924 w = velRadial[0,:]
4716 4925 winds = velRadial.copy()
4717 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 4932 if my_nbeams == 2:
4933
4720 4934 zenith = numpy.array(zenith)
4721 4935 zenith -= zenithCorrection
4722 4936 zenith *= numpy.pi/180
@@ -4734,25 +4948,73 class EWDriftsEstimation(Operation):
4734 4948 w_e = velRadial1[1,:]
4735 4949 w_w_err = velRerr[0,:]
4736 4950 w_e_err = velRerr[1,:]
4737
4738 val = (numpy.isfinite(w_w)==False).nonzero()
4739 val = val[0]
4740 bad = val
4741 if len(bad) > 0 :
4742 w_w[bad] = w_wC[bad]
4743 w_w_err[bad]= numpy.nan
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,:]
4955
4956 tini=time.localtime(dataOut.utctime)
4957 #print(tini[3],tini[4])
4958 #val = (numpy.isfinite(w_w)==False).nonzero()
4959 #val = val[0]
4960 #bad = val
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 4976 smooth_eC=ebufc[4,:]
4745 4977 p_e0C = rbufc[4,:]
4746 4978 p_e1C = rbufc[5,:]
4747 4979 w_eC = rbufc[6,:]*-1
4748 4980 t_eC = rbufc[7,:]
4749 val = (numpy.isfinite(w_e)==False).nonzero()
4750 val = val[0]
4751 bad = val
4752 if len(bad) > 0 :
4753 w_e[bad] = w_eC[bad]
4754 w_e_err[bad]= numpy.nan
4755
4981
4982 val = (numpy.isfinite(p_e0)==False).nonzero()
4983 p_e0[val]=0
4984 val = (numpy.isfinite(p_e1)==False).nonzero()
4985 p_e1[val]=0
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
5017
4756 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 5019 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
4758 5020
@@ -4760,13 +5022,22 class EWDriftsEstimation(Operation):
4760 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 5024 winds = numpy.vstack((w,u))
4763
5025 #winds = numpy.vstack((w,u,w_err,u_err))
4764 5026 dataOut.heightList = heiRang1
4765 snr1 = 10*numpy.log10(SNR1[0])
5027 #snr1 = 10*numpy.log10(SNR1[0])
4766 5028 dataOut.data_output = winds
4767 #snr1 = 10*numpy.log10(SNR1[0])# estaba comentado
4768 dataOut.data_snr1 = numpy.reshape(snr1,(1,snr1.shape[0]))
4769 #print("data_snr1",dataOut.data_snr1)
5029 range1 = dataOut.heightList
5030 nhei = len(range1)
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 5041 dataOut.utctimeInit = dataOut.utctime
4771 5042 dataOut.outputInterval = dataOut.timeInterval
4772 5043
@@ -4774,8 +5045,6 class EWDriftsEstimation(Operation):
4774 5045 jrange = 450 #900 para HA drifts
4775 5046 deltah = 15.0 #dataOut.spacing(0) 25 HAD
4776 5047 h0 = 0.0 #dataOut.first_height(0)
4777 heights = dataOut.heightList
4778 nhei = len(heights)
4779 5048
4780 5049 range1 = numpy.arange(nhei) * deltah + h0
4781 5050 jhei = (range1 >= hei_aver0).nonzero()
@@ -4840,18 +5109,95 class EWDriftsEstimation(Operation):
4840 5109 uA = uA[6*h_avgs:nhei_avg-0]
4841 5110 uA_err = uA_err[6*h_avgs:nhei_avg-0]
4842 5111 dataOut.drifts_avg = numpy.vstack((wA,uA))
5112
4843 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 5120 tini=time.localtime(dataOut.utctime)
4845 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 5125 f1 = open(nfile,'a')
4848 5126 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
4849 5127 driftavgstr=str(dataOut.drifts_avg)
4850 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')
5132 f1.close()
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')
4852 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):
4853 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)
4854 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)
5198 return dataOut
5199
5200
4855 5201
4856 5202 #--------------- Non Specular Meteor ----------------
4857 5203
General Comments 0
You need to be logged in to leave comments. Login now