##// END OF EJS Templates
schain
RSS
LAST UPDATE AFTER PERCY & ALEX REVISION
Alexander Valdez -
r1672:65c1739783d6
parent child
Show More
Side by Side Unified
Context file:
r1672:65c1739783d6 Loading diff...:
Show file before | Show file after | Hide comments
@@ -85,10 +85,11 class ParametersProc(ProcessingUnit):
85 self.dataOut.timeInterval1 = self.dataIn.timeInterval
85 self.dataOut.timeInterval1 = self.dataIn.timeInterval
86 self.dataOut.heightList = self.dataIn.heightList
86 self.dataOut.heightList = self.dataIn.heightList
87 self.dataOut.frequency = self.dataIn.frequency
87 self.dataOut.frequency = self.dataIn.frequency
88 self.dataOut.runNextUnit = self.dataIn.runNextUnit
88
89
90 def run(self, runNextUnit=0):
89
91
90 def run(self):
92 self.dataIn.runNextUnit = runNextUnit
91
92 #---------------------- Voltage Data ---------------------------
93 #---------------------- Voltage Data ---------------------------
93
94
94 if self.dataIn.type == "Voltage":
95 if self.dataIn.type == "Voltage":
@@ -1454,7 +1455,6 class SpectralMoments(Operation):
1454 dataOut.data_snr = data_param[:,0]
1455 dataOut.data_snr = data_param[:,0]
1455 dataOut.data_pow = data_param[:,6] # to compare with type0 proccessing
1456 dataOut.data_pow = data_param[:,6] # to compare with type0 proccessing
1456 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, data_param[:,3], data_param[:,4],data_param[:,5]),axis=2)
1457 dataOut.spcpar=numpy.stack((dataOut.data_dop,dataOut.data_width,dataOut.data_snr, data_param[:,3], data_param[:,4],data_param[:,5]),axis=2)
1457
1458 else:
1458 else:
1459 dataOut.moments = data_param[:,1:,:]
1459 dataOut.moments = data_param[:,1:,:]
1460 dataOut.data_snr = data_param[:,0]
1460 dataOut.data_snr = data_param[:,0]
@@ -2223,6 +2223,7 class SpectralFitting(Operation):
2223 return moments
2223 return moments
2224
2224
2225 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
2225 def __DiffCoherent(self, spectra, cspectra, dataOut, noise, snrth, coh_th, hei_th):
2226
2226 nProf = dataOut.nProfiles
2227 nProf = dataOut.nProfiles
2227 heights = dataOut.heightList
2228 heights = dataOut.heightList
2228 nHei = len(heights)
2229 nHei = len(heights)
@@ -2248,20 +2249,17 class SpectralFitting(Operation):
2248
2249
2249 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
2250 if coh_th == None : coh_th = numpy.array([0.75,0.65,0.15]) # 0.65
2250 if hei_th == None : hei_th = numpy.array([60,300,650])
2251 if hei_th == None : hei_th = numpy.array([60,300,650])
2251 for ic in range(2):
2252 for ic in range(nPairs):
2252 pair = crosspairs[ic]
2253 pair = crosspairs[ic]
2253 #si el SNR es mayor que el SNR threshold los datos se toman coherentes
2254 #si el SNR es mayor que el SNR threshold los datos se toman coherentes
2254 s_n0 = power[pair[0],:]/noise[pair[0]]
2255 s_n0 = power[pair[0],:]/noise[pair[0]]
2255 s_n1 = power[pair[1],:]/noise[pair[1]]
2256 s_n1 = power[pair[1],:]/noise[pair[1]]
2256
2257 valid1 =(s_n0>=snr_th).nonzero()
2257 valid1 =(s_n0>=snr_th).nonzero()
2258 valid2 = (s_n1>=snr_th).nonzero()
2258 valid2 = (s_n1>=snr_th).nonzero()
2259 #valid = valid2 + valid1 #numpy.concatenate((valid1,valid2), axis=None)
2260 valid1 = numpy.array(valid1[0])
2259 valid1 = numpy.array(valid1[0])
2261 valid2 = numpy.array(valid2[0])
2260 valid2 = numpy.array(valid2[0])
2262 valid = valid1
2261 valid = valid1
2263 for iv in range(len(valid2)):
2262 for iv in range(len(valid2)):
2264 #for ivv in range(len(valid1)) :
2265 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2263 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2266 if len(indv[0]) == 0 :
2264 if len(indv[0]) == 0 :
2267 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2265 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
@@ -2269,17 +2267,13 class SpectralFitting(Operation):
2269 my_coh_aver[pair[0],valid]=1
2267 my_coh_aver[pair[0],valid]=1
2270 my_coh_aver[pair[1],valid]=1
2268 my_coh_aver[pair[1],valid]=1
2271 # si la coherencia es mayor a la coherencia threshold los datos se toman
2269 # si la coherencia es mayor a la coherencia threshold los datos se toman
2272 #print my_coh_aver[0,:]
2273 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)))
2270 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)))
2274 #print('coh',numpy.absolute(coh))
2275 for ih in range(len(hei_th)):
2271 for ih in range(len(hei_th)):
2276 hvalid = (heights>hei_th[ih]).nonzero()
2272 hvalid = (heights>hei_th[ih]).nonzero()
2277 hvalid = hvalid[0]
2273 hvalid = hvalid[0]
2278 if len(hvalid)>0:
2274 if len(hvalid)>0:
2279 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
2275 valid = (numpy.absolute(coh[hvalid])>coh_th[ih]).nonzero()
2280 valid = valid[0]
2276 valid = valid[0]
2281 #print('hvalid:',hvalid)
2282 #print('valid', valid)
2283 if len(valid)>0:
2277 if len(valid)>0:
2284 my_coh_aver[pair[0],hvalid[valid]] =1
2278 my_coh_aver[pair[0],hvalid[valid]] =1
2285 my_coh_aver[pair[1],hvalid[valid]] =1
2279 my_coh_aver[pair[1],hvalid[valid]] =1
@@ -2295,7 +2289,7 class SpectralFitting(Operation):
2295 my_incoh_aver[pair[1],incoh_echoes] = 1
2289 my_incoh_aver[pair[1],incoh_echoes] = 1
2296
2290
2297
2291
2298 for ic in range(2):
2292 for ic in range(nPairs):
2299 pair = crosspairs[ic]
2293 pair = crosspairs[ic]
2300
2294
2301 valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
2295 valid1 =(my_coh_aver[pair[0],:]==1 ).nonzero()
@@ -2303,29 +2297,24 class SpectralFitting(Operation):
2303 valid1 = numpy.array(valid1[0])
2297 valid1 = numpy.array(valid1[0])
2304 valid2 = numpy.array(valid2[0])
2298 valid2 = numpy.array(valid2[0])
2305 valid = valid1
2299 valid = valid1
2306 #print valid1 , valid2
2300
2307 for iv in range(len(valid2)):
2301 for iv in range(len(valid2)):
2308 #for ivv in range(len(valid1)) :
2302
2309 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2303 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2310 if len(indv[0]) == 0 :
2304 if len(indv[0]) == 0 :
2311 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2305 valid = numpy.concatenate((valid,valid2[iv]), axis=None)
2312 #print valid
2313 #valid = numpy.concatenate((valid1,valid2), axis=None)
2314 valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
2306 valid1 =(my_coh_aver[pair[0],:] !=1 ).nonzero()
2315 valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
2307 valid2 = (my_coh_aver[pair[1],:] !=1).nonzero()
2316 valid1 = numpy.array(valid1[0])
2308 valid1 = numpy.array(valid1[0])
2317 valid2 = numpy.array(valid2[0])
2309 valid2 = numpy.array(valid2[0])
2318 incoh_echoes = valid1
2310 incoh_echoes = valid1
2319 #print valid1, valid2
2320 #incoh_echoes= numpy.concatenate((valid1,valid2), axis=None)
2321 for iv in range(len(valid2)):
2311 for iv in range(len(valid2)):
2322 #for ivv in range(len(valid1)) :
2312
2323 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2313 indv = numpy.array((valid1 == valid2[iv]).nonzero())
2324 if len(indv[0]) == 0 :
2314 if len(indv[0]) == 0 :
2325 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
2315 incoh_echoes = numpy.concatenate(( incoh_echoes,valid2[iv]), axis=None)
2326 #print incoh_echoes
2316
2327 if len(valid)>0:
2317 if len(valid)>0:
2328 #print pair
2329 coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
2318 coh_spectra[pair[0],:,valid] = spectra[pair[0],:,valid]
2330 coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
2319 coh_spectra[pair[1],:,valid] = spectra[pair[1],:,valid]
2331 coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
2320 coh_cspectra[ic,:,valid] = cspectra[ic,:,valid]
@@ -2339,6 +2328,7 class SpectralFitting(Operation):
2339 incoh_aver[pair[1],incoh_echoes]=1
2328 incoh_aver[pair[1],incoh_echoes]=1
2340 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
2329 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
2341
2330
2331
2342 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
2332 def __CleanCoherent(self,snrth, spectra, cspectra, coh_aver,dataOut, noise,clean_coh_echoes,index):
2343
2333
2344 nProf = dataOut.nProfiles
2334 nProf = dataOut.nProfiles
@@ -2349,10 +2339,7 class SpectralFitting(Operation):
2349 crosspairs = dataOut.groupList
2339 crosspairs = dataOut.groupList
2350 nPairs = len(crosspairs)
2340 nPairs = len(crosspairs)
2351
2341
2352 #data = dataOut.data_pre[0]
2353 absc = dataOut.abscissaList[:-1]
2342 absc = dataOut.abscissaList[:-1]
2354 #noise = dataOut.noise
2355 #nChannel = data.shape[0]
2356 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
2343 data_param = numpy.zeros((nChan, 4, spectra.shape[2]))
2357 clean_coh_spectra = spectra.copy()
2344 clean_coh_spectra = spectra.copy()
2358 clean_coh_cspectra = cspectra.copy()
2345 clean_coh_cspectra = cspectra.copy()
@@ -2364,17 +2351,12 class SpectralFitting(Operation):
2364 rtime0 = [6,18] # periodo sin ESF
2351 rtime0 = [6,18] # periodo sin ESF
2365 rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
2352 rtime1 = [10.5,13.5] # periodo con alta coherencia y alto ancho espectral (esperado): SOL.
2366
2353
2367 time = index*5./60
2354 time = index*5./60 # en base a 5 min de proceso
2368 if clean_coh_echoes == 1 :
2355 if clean_coh_echoes == 1 :
2369 for ind in range(nChan):
2356 for ind in range(nChan):
2370 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
2357 data_param[ind,:,:] = self.__calculateMoments( spectra[ind,:,:] , absc , noise[ind] )
2371 #print data_param[:,3]
2372 spwd = data_param[:,3]
2358 spwd = data_param[:,3]
2373 #print spwd.shape
2374 # 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
2359 # 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
2375 #spwd1=[ 1.65607, 1.43416, 0.500373, 0.208361, 0.000000, 26.7767, 22.5936, 26.7530, 20.6962, 29.1098, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 28.0300, 27.0511, 27.8810, 26.3126, 27.8445, 24.6181, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000]
2376 #spwd=numpy.array([spwd1,spwd1,spwd1,spwd1])
2377 #print spwd.shape, heights.shape,coh_aver.shape
2378 # para obtener spwd
2360 # para obtener spwd
2379 for ic in range(nPairs):
2361 for ic in range(nPairs):
2380 pair = crosspairs[ic]
2362 pair = crosspairs[ic]
@@ -2399,23 +2381,15 class SpectralFitting(Operation):
2399
2381
2400 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
2382 return clean_coh_spectra, clean_coh_cspectra, clean_coh_aver
2401
2383
2402 #def CleanRayleigh(self,dataOut,spectra,cspectra,out_spectra,out_cspectra,sat_spectra,sat_cspectra,crosspairs,heights, channels, nProf,nHei,nChan,nPairs,nIncohInt,nBlocks):
2403 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
2384 def CleanRayleigh(self,dataOut,spectra,cspectra,save_drifts):
2404 #import matplotlib.pyplot as plt
2385
2405 #for k in range(149):
2386 rfunc = cspectra.copy()
2406
2407 # self.bloque0[:,:,:,k] = spectra[:,:,0:nHei]
2408 # self.bloques[:,:,:,k] = cspectra[:,:,0:nHei]
2409 #if self.i==nBlocks:
2410 # self.i==0
2411 rfunc = cspectra.copy() #self.bloques
2412 n_funct = len(rfunc[0,:,0,0])
2387 n_funct = len(rfunc[0,:,0,0])
2413 val_spc = spectra*0.0 #self.bloque0*0.0
2388 val_spc = spectra*0.0
2414 val_cspc = cspectra*0.0 #self.bloques*0.0
2389 val_cspc = cspectra*0.0
2415 in_sat_spectra = spectra.copy() #self.bloque0
2390 in_sat_spectra = spectra.copy()
2416 in_sat_cspectra = cspectra.copy() #self.bloques
2391 in_sat_cspectra = cspectra.copy()
2417
2392
2418 #print( rfunc.shape)
2419 min_hei = 200
2393 min_hei = 200
2420 nProf = dataOut.nProfiles
2394 nProf = dataOut.nProfiles
2421 heights = dataOut.heightList
2395 heights = dataOut.heightList
@@ -2426,20 +2400,17 class SpectralFitting(Operation):
2426 nPairs = len(crosspairs)
2400 nPairs = len(crosspairs)
2427 hval=(heights >= min_hei).nonzero()
2401 hval=(heights >= min_hei).nonzero()
2428 ih=hval[0]
2402 ih=hval[0]
2429 #print numpy.absolute(rfunc[:,0,0,14])
2430 for ih in range(hval[0][0],nHei):
2403 for ih in range(hval[0][0],nHei):
2431 for ifreq in range(nProf):
2404 for ifreq in range(nProf):
2432 for ii in range(n_funct):
2405 for ii in range(n_funct):
2433
2406
2434 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
2407 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))
2435 #print numpy.amin(func2clean)
2436 val = (numpy.isfinite(func2clean)==True).nonzero()
2408 val = (numpy.isfinite(func2clean)==True).nonzero()
2437 if len(val)>0:
2409 if len(val)>0:
2438 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
2410 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
2439 if min_val <= -40 : min_val = -40
2411 if min_val <= -40 : min_val = -40
2440 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
2412 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
2441 if max_val >= 200 : max_val = 200
2413 if max_val >= 200 : max_val = 200
2442 #print min_val, max_val
2443 step = 1
2414 step = 1
2444 #Getting bins and the histogram
2415 #Getting bins and the histogram
2445 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
2416 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
@@ -2454,9 +2425,6 class SpectralFitting(Operation):
2454 except:
2425 except:
2455 mode = mean
2426 mode = mean
2456 stdv = sigma
2427 stdv = sigma
2457 # 7.84616 53.9307 3.61863
2458 #stdv = 3.61863 # 2.99089
2459 #mode = 53.9307 #7.79008
2460
2428
2461 #Removing echoes greater than mode + 3*stdv
2429 #Removing echoes greater than mode + 3*stdv
2462 factor_stdv = 2.5
2430 factor_stdv = 2.5
@@ -2467,38 +2435,14 class SpectralFitting(Operation):
2467 cross_pairs = crosspairs[ii]
2435 cross_pairs = crosspairs[ii]
2468 #Getting coherent echoes which are removed.
2436 #Getting coherent echoes which are removed.
2469 if len(novall[0]) > 0:
2437 if len(novall[0]) > 0:
2470 #val_spc[(0,1),novall[a],ih] = 1
2471 #val_spc[,(2,3),novall[a],ih] = 1
2472 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
2438 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
2473 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
2439 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
2474 val_cspc[novall[0],ii,ifreq,ih] = 1
2440 val_cspc[novall[0],ii,ifreq,ih] = 1
2475 #print("OUT NOVALL 1")
2476 #Removing coherent from ISR data
2441 #Removing coherent from ISR data
2477 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
2442 spectra[noval,cross_pairs[0],ifreq,ih] = numpy.nan
2478 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
2443 spectra[noval,cross_pairs[1],ifreq,ih] = numpy.nan
2479 cspectra[noval,ii,ifreq,ih] = numpy.nan
2444 cspectra[noval,ii,ifreq,ih] = numpy.nan
2480 #no sale es para savedrifts >2
2481 ''' channels = channels
2482 cross_pairs = cross_pairs
2483 #print("OUT NOVALL 2")
2484
2445
2485 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
2486 vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
2487 vcross = numpy.concatenate((vcross0,vcross1),axis=None)
2488 #print('vcros =', vcross)
2489
2490 #Getting coherent echoes which are removed.
2491 if len(novall) > 0:
2492 #val_spc[novall,ii,ifreq,ih] = 1
2493 val_spc[ii,ifreq,ih,novall] = 1
2494 if len(vcross) > 0:
2495 val_cspc[vcross,ifreq,ih,novall] = 1
2496
2497 #Removing coherent from ISR data.
2498 self.bloque0[ii,ifreq,ih,noval] = numpy.nan
2499 if len(vcross) > 0:
2500 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
2501 '''
2502 #Getting average of the spectra and cross-spectra from incoherent echoes.
2446 #Getting average of the spectra and cross-spectra from incoherent echoes.
2503 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
2447 out_spectra = numpy.zeros([nChan,nProf,nHei], dtype=float) #+numpy.nan
2504 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
2448 out_cspectra = numpy.zeros([nPairs,nProf,nHei], dtype=complex) #+numpy.nan
@@ -2507,18 +2451,14 class SpectralFitting(Operation):
2507 for ich in range(nChan):
2451 for ich in range(nChan):
2508 tmp = spectra[:,ich,ifreq,ih]
2452 tmp = spectra[:,ich,ifreq,ih]
2509 valid = (numpy.isfinite(tmp[:])==True).nonzero()
2453 valid = (numpy.isfinite(tmp[:])==True).nonzero()
2510 #print('TMP',tmp)
2511 if len(valid[0]) >0 :
2454 if len(valid[0]) >0 :
2512 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2455 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2513 #for icr in range(nPairs):
2514 for icr in range(nPairs):
2456 for icr in range(nPairs):
2515 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
2457 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
2516 valid = (numpy.isfinite(tmp)==True).nonzero()
2458 valid = (numpy.isfinite(tmp)==True).nonzero()
2517 if len(valid[0]) > 0:
2459 if len(valid[0]) > 0:
2518 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2460 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2519 # print('##########################################################')
2520 #Removing fake coherent echoes (at least 4 points around the point)
2461 #Removing fake coherent echoes (at least 4 points around the point)
2521
2522 val_spectra = numpy.sum(val_spc,0)
2462 val_spectra = numpy.sum(val_spc,0)
2523 val_cspectra = numpy.sum(val_cspc,0)
2463 val_cspectra = numpy.sum(val_cspc,0)
2524
2464
@@ -2564,39 +2504,26 class SpectralFitting(Operation):
2564 valid = (numpy.isfinite(tmp)).nonzero()
2504 valid = (numpy.isfinite(tmp)).nonzero()
2565 if len(valid[0]) > 0:
2505 if len(valid[0]) > 0:
2566 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2506 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
2567 #self.__dataReady= True
2568 #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
2569 #if not self.__dataReady:
2570 #return None, None
2571 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
2507 return out_spectra, out_cspectra,sat_spectra,sat_cspectra
2572 def REM_ISOLATED_POINTS(self,array,rth):
2508 def REM_ISOLATED_POINTS(self,array,rth):
2573 if rth == None : rth = 4
2509 if rth == None : rth = 4
2574
2575 num_prof = len(array[0,:,0])
2510 num_prof = len(array[0,:,0])
2576 num_hei = len(array[0,0,:])
2511 num_hei = len(array[0,0,:])
2577 n2d = len(array[:,0,0])
2512 n2d = len(array[:,0,0])
2578
2513
2579 for ii in range(n2d) :
2514 for ii in range(n2d) :
2580 #print ii,n2d
2581 tmp = array[ii,:,:]
2515 tmp = array[ii,:,:]
2582 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
2583 tmp = numpy.reshape(tmp,num_prof*num_hei)
2516 tmp = numpy.reshape(tmp,num_prof*num_hei)
2584 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
2517 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
2585 indxs2 = (tmp > 0).nonzero()
2518 indxs2 = (tmp > 0).nonzero()
2586
2587 indxs1 = (indxs1[0])
2519 indxs1 = (indxs1[0])
2588 indxs2 = indxs2[0]
2520 indxs2 = indxs2[0]
2589 #indxs1 = numpy.array(indxs1[0])
2590 #indxs2 = numpy.array(indxs2[0])
2591 indxs = None
2521 indxs = None
2592 #print indxs1 , indxs2
2593 for iv in range(len(indxs2)):
2522 for iv in range(len(indxs2)):
2594 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
2523 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
2595 #print len(indxs2), indv
2596 if len(indv[0]) > 0 :
2524 if len(indv[0]) > 0 :
2597 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
2525 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
2598 indxs = indxs[1:]
2526 indxs = indxs[1:]
2599 #print indxs, len(indxs)
2600 if len(indxs) < 4 :
2527 if len(indxs) < 4 :
2601 array[ii,:,:] = 0.
2528 array[ii,:,:] = 0.
2602 return
2529 return
@@ -2604,41 +2531,30 class SpectralFitting(Operation):
2604 xpos = numpy.mod(indxs ,num_hei)
2531 xpos = numpy.mod(indxs ,num_hei)
2605 ypos = (indxs / num_hei)
2532 ypos = (indxs / num_hei)
2606 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
2533 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
2607 #print sx
2608 xpos = xpos[sx]
2534 xpos = xpos[sx]
2609 ypos = ypos[sx]
2535 ypos = ypos[sx]
2610 # *********************************** Cleaning isolated points **********************************
2536 # *********************************** Cleaning isolated points **********************************
2611 ic = 0
2537 ic = 0
2612 while True :
2538 while True :
2613 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
2539 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
2614 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
2615 #plt.plot(r)
2616 #plt.show()
2617 no_coh1 = (numpy.isfinite(r)==True).nonzero()
2540 no_coh1 = (numpy.isfinite(r)==True).nonzero()
2618 no_coh2 = (r <= rth).nonzero()
2541 no_coh2 = (r <= rth).nonzero()
2619 #print r, no_coh1, no_coh2
2620 no_coh1 = numpy.array(no_coh1[0])
2542 no_coh1 = numpy.array(no_coh1[0])
2621 no_coh2 = numpy.array(no_coh2[0])
2543 no_coh2 = numpy.array(no_coh2[0])
2622 no_coh = None
2544 no_coh = None
2623 #print valid1 , valid2
2624 for iv in range(len(no_coh2)):
2545 for iv in range(len(no_coh2)):
2625 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
2546 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
2626 if len(indv[0]) > 0 :
2547 if len(indv[0]) > 0 :
2627 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
2548 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
2628 no_coh = no_coh[1:]
2549 no_coh = no_coh[1:]
2629 #print len(no_coh), no_coh
2630 if len(no_coh) < 4 :
2550 if len(no_coh) < 4 :
2631 #print xpos[ic], ypos[ic], ic
2632 xpos[ic] = numpy.nan
2551 xpos[ic] = numpy.nan
2633 ypos[ic] = numpy.nan
2552 ypos[ic] = numpy.nan
2634
2553
2635 ic = ic + 1
2554 ic = ic + 1
2636 if (ic == len(indxs)) :
2555 if (ic == len(indxs)) :
2637 break
2556 break
2638 #print( xpos, ypos)
2639
2640 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
2557 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
2641 #print indxs[0]
2642 if len(indxs[0]) < 4 :
2558 if len(indxs[0]) < 4 :
2643 array[ii,:,:] = 0.
2559 array[ii,:,:] = 0.
2644 return
2560 return
@@ -2652,20 +2568,12 class SpectralFitting(Operation):
2652
2568
2653 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
2569 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
2654 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
2570 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
2655
2656 #print array.shape
2657 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
2658 #print tmp.shape
2659
2660 return array
2571 return array
2572
2661 def moments(self,doppler,yarray,npoints):
2573 def moments(self,doppler,yarray,npoints):
2662 ytemp = yarray
2574 ytemp = yarray
2663 #val = WHERE(ytemp GT 0,cval)
2664 #if cval == 0 : val = range(npoints-1)
2665 val = (ytemp > 0).nonzero()
2575 val = (ytemp > 0).nonzero()
2666 val = val[0]
2576 val = val[0]
2667 #print('hvalid:',hvalid)
2668 #print('valid', valid)
2669 if len(val) == 0 : val = range(npoints-1)
2577 if len(val) == 0 : val = range(npoints-1)
2670
2578
2671 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
2579 ynew = 0.5*(ytemp[val[0]]+ytemp[val[len(val)-1]])
@@ -2684,15 +2592,15 class SpectralFitting(Operation):
2684
2592
2685
2593
2686
2594
2687
2595 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):
2688
2596 if not numpy.any(proc):
2689
2690 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None, filec=None,coh_th=None, hei_th=None):
2691 nChannels = dataOut.nChannels
2597 nChannels = dataOut.nChannels
2692 nHeights= dataOut.heightList.size
2598 nHeights= dataOut.heightList.size
2693 nProf = dataOut.nProfiles
2599 nProf = dataOut.nProfiles
2600 if numpy.any(taver): taver=int(taver)
2601 else : taver = 5
2694 tini=time.localtime(dataOut.utctime)
2602 tini=time.localtime(dataOut.utctime)
2695 if (tini.tm_min % 5) == 0 and (tini.tm_sec < 5 and self.fint==0):
2603 if (tini.tm_min % taver) == 0 and (tini.tm_sec < 5 and self.fint==0):
2696 self.index = 0
2604 self.index = 0
2697 jspc = self.buffer
2605 jspc = self.buffer
2698 jcspc = self.buffer2
2606 jcspc = self.buffer2
@@ -2702,14 +2610,14 class SpectralFitting(Operation):
2702 self.buffer3 = dataOut.noise
2610 self.buffer3 = dataOut.noise
2703 self.fint = 1
2611 self.fint = 1
2704 if numpy.any(jspc) :
2612 if numpy.any(jspc) :
2705 jspc= numpy.reshape(jspc,(int(len(jspc)/4),nChannels,nProf,nHeights))
2613 jspc= numpy.reshape(jspc,(int(len(jspc)/nChannels),nChannels,nProf,nHeights))
2706 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/2),2,nProf,nHeights))
2614 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/int(nChannels/2)),int(nChannels/2),nProf,nHeights))
2707 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/4),nChannels))
2615 jnoise= numpy.reshape(jnoise,(int(len(jnoise)/nChannels),nChannels))
2708 else:
2616 else:
2709 dataOut.flagNoData = True
2617 dataOut.flagNoData = True
2710 return dataOut
2618 return dataOut
2711 else :
2619 else :
2712 if (tini.tm_min % 5) == 0 : self.fint = 1
2620 if (tini.tm_min % taver) == 0 : self.fint = 1
2713 else : self.fint = 0
2621 else : self.fint = 0
2714 self.index += 1
2622 self.index += 1
2715 if numpy.any(self.buffer):
2623 if numpy.any(self.buffer):
@@ -2725,12 +2633,13 class SpectralFitting(Operation):
2725 if path != None:
2633 if path != None:
2726 sys.path.append(path)
2634 sys.path.append(path)
2727 self.library = importlib.import_module(file)
2635 self.library = importlib.import_module(file)
2636 if filec != None:
2637 self.weightf = importlib.import_module(filec)
2728
2638
2729 #To be inserted as a parameter
2639 #To be inserted as a parameter
2730 groupArray = numpy.array(groupList)
2640 groupArray = numpy.array(groupList)
2731 #groupArray = numpy.array([[0,1],[2,3]])
2641 #groupArray = numpy.array([[0,1],[2,3]])
2732 dataOut.groupList = groupArray
2642 dataOut.groupList = groupArray
2733
2734 nGroups = groupArray.shape[0]
2643 nGroups = groupArray.shape[0]
2735 nChannels = dataOut.nChannels
2644 nChannels = dataOut.nChannels
2736 nHeights = dataOut.heightList.size
2645 nHeights = dataOut.heightList.size
@@ -2738,6 +2647,13 class SpectralFitting(Operation):
2738 #Parameters Array
2647 #Parameters Array
2739 dataOut.data_param = None
2648 dataOut.data_param = None
2740 dataOut.data_paramC = None
2649 dataOut.data_paramC = None
2650 dataOut.clean_num_aver = None
2651 dataOut.coh_num_aver = None
2652 dataOut.tmp_spectra_i = None
2653 dataOut.tmp_cspectra_i = None
2654 dataOut.tmp_spectra_c = None
2655 dataOut.tmp_cspectra_c = None
2656 dataOut.index = None
2741
2657
2742 #Set constants
2658 #Set constants
2743 constants = self.library.setConstants(dataOut)
2659 constants = self.library.setConstants(dataOut)
@@ -2747,7 +2663,6 class SpectralFitting(Operation):
2747 ippSeconds = dataOut.ippSeconds
2663 ippSeconds = dataOut.ippSeconds
2748 K = dataOut.nIncohInt
2664 K = dataOut.nIncohInt
2749 pairsArray = numpy.array(dataOut.pairsList)
2665 pairsArray = numpy.array(dataOut.pairsList)
2750
2751 snrth= 20
2666 snrth= 20
2752 spectra = dataOut.data_spc
2667 spectra = dataOut.data_spc
2753 cspectra = dataOut.data_cspc
2668 cspectra = dataOut.data_cspc
@@ -2768,7 +2683,8 class SpectralFitting(Operation):
2768 if not self.isConfig:
2683 if not self.isConfig:
2769 self.isConfig = True
2684 self.isConfig = True
2770
2685
2771 index = tini.tm_hour*12+tini.tm_min/5
2686 index = tini.tm_hour*12+tini.tm_min/taver
2687 dataOut.index= index
2772 jspc = jspc/N/N
2688 jspc = jspc/N/N
2773 jcspc = jcspc/N/N
2689 jcspc = jcspc/N/N
2774 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
2690 tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.CleanRayleigh(dataOut,jspc,jcspc,2)
@@ -2778,17 +2694,68 class SpectralFitting(Operation):
2778 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
2694 clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
2779 dataOut.data_spc = incoh_spectra
2695 dataOut.data_spc = incoh_spectra
2780 dataOut.data_cspc = incoh_cspectra
2696 dataOut.data_cspc = incoh_cspectra
2781
2782 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
2697 clean_num_aver = incoh_aver*len(jspc[:,0,0,0])
2783 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
2698 coh_num_aver = clean_coh_aver*len(jspc[:,0,0,0])
2699 dataOut.clean_num_aver = clean_num_aver
2700 dataOut.coh_num_aver = coh_num_aver
2701 dataOut.tmp_spectra_i = incoh_spectra
2702 dataOut.tmp_cspectra_i = incoh_cspectra
2703 dataOut.tmp_spectra_c = clean_coh_spectra
2704 dataOut.tmp_cspectra_c = clean_coh_cspectra
2784 #List of possible combinations
2705 #List of possible combinations
2785 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
2706 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
2786 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
2707 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
2787
2788 if getSNR:
2708 if getSNR:
2789 listChannels = groupArray.reshape((groupArray.size))
2709 listChannels = groupArray.reshape((groupArray.size))
2790 listChannels.sort()
2710 listChannels.sort()
2791 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
2711 dataOut.data_SNR = self.__getSNR(dataOut.data_spc[listChannels,:,:], noise[listChannels])
2712 else:
2713 clean_num_aver = dataOut.clean_num_aver
2714 coh_num_aver = dataOut.coh_num_aver
2715 dataOut.data_spc = dataOut.tmp_spectra_i
2716 dataOut.data_cspc = dataOut.tmp_cspectra_i
2717 clean_coh_spectra = dataOut.tmp_spectra_c
2718 clean_coh_cspectra = dataOut.tmp_cspectra_c
2719 jspectra = dataOut.data_spc+clean_coh_spectra
2720 nHeights = len(dataOut.heightList) # nhei
2721 nProf = int(dataOut.nProfiles)
2722 dataOut.nProfiles = nProf
2723 dataOut.data_param = None
2724 dataOut.data_paramC = None
2725 dataOut.code = numpy.array([[-1.,-1.,1.],[1.,1.,-1.]])
2726 #M=600
2727 #N=200
2728 dataOut.flagDecodeData=True
2729 M = int(dataOut.normFactor)
2730 N = int(dataOut.nFFTPoints)
2731 dataOut.nFFTPoints = N
2732 dataOut.nIncohInt= int(dataOut.nIncohInt)
2733 dataOut.nProfiles = int(dataOut.nProfiles)
2734 dataOut.nCohInt = int(dataOut.nCohInt)
2735 print('sale',dataOut.nProfiles,dataOut.nHeights)
2736 #dataOut.nFFTPoints=nprofs
2737 #dataOut.normFactor = nprofs
2738 dataOut.channelList = channelList
2739 #dataOut.ippFactor=1
2740 #ipp = ipp/150*1.e-3
2741 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
2742 #dataOut.ippSeconds=ipp
2743 absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
2744 print('sale 2',dataOut.ippSeconds,M,N)
2745 print('Empieza procesamiento offline')
2746 if path != None:
2747 sys.path.append(path)
2748 self.library = importlib.import_module(file)
2749 constants = self.library.setConstants(dataOut)
2750 constants['M'] = M
2751 dataOut.constants = constants
2752
2753 groupArray = numpy.array(groupList)
2754 dataOut.groupList = groupArray
2755 nGroups = groupArray.shape[0]
2756 #List of possible combinations
2757 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
2758 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
2792 if dataOut.data_paramC is None:
2759 if dataOut.data_paramC is None:
2793 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
2760 dataOut.data_paramC = numpy.zeros((nGroups*4, nHeights,2))*numpy.nan
2794 for i in range(nGroups):
2761 for i in range(nGroups):
@@ -2807,17 +2774,18 class SpectralFitting(Operation):
2807 dataCross = dataCross**2
2774 dataCross = dataCross**2
2808 nhei = nHeights
2775 nhei = nHeights
2809 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
2776 poweri = numpy.sum(dataOut.data_spc[:,1:nProf-0,:],axis=1)/clean_num_aver[:,:]
2810 if i == 0 : my_noises = numpy.zeros(4,dtype=float) #FLTARR(4)
2777 if i == 0 : my_noises = numpy.zeros(4,dtype=float)
2811 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
2778 n0i = numpy.nanmin(poweri[0+i*2,0:nhei-0])/(nProf-1)
2812 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
2779 n1i = numpy.nanmin(poweri[1+i*2,0:nhei-0])/(nProf-1)
2813 n0 = n0i
2780 n0 = n0i
2814 n1= n1i
2781 n1= n1i
2815 my_noises[2*i+0] = n0
2782 my_noises[2*i+0] = n0
2816 my_noises[2*i+1] = n1
2783 my_noises[2*i+1] = n1
2817 snrth = -15.0 #-16 -10
2784 snrth = -15.0 # -4 -16 -25
2818 snrth = 10**(snrth/10.0)
2785 snrth = 10**(snrth/10.0)
2819 jvelr = numpy.zeros(nHeights, dtype = 'float')
2786 jvelr = numpy.zeros(nHeights, dtype = 'float')
2820 hvalid = [0]
2787 hvalid = [0]
2788 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,:])
2821 for h in range(nHeights):
2789 for h in range(nHeights):
2822 smooth = clean_num_aver[i+1,h]
2790 smooth = clean_num_aver[i+1,h]
2823 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
2791 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
@@ -2826,15 +2794,23 class SpectralFitting(Operation):
2826 signal1 = signalpn1-n1
2794 signal1 = signalpn1-n1
2827 snr0 = numpy.sum(signal0/n0)/(nProf-1)
2795 snr0 = numpy.sum(signal0/n0)/(nProf-1)
2828 snr1 = numpy.sum(signal1/n1)/(nProf-1)
2796 snr1 = numpy.sum(signal1/n1)/(nProf-1)
2797 gamma = coh2[:,h]
2798 indxs = (numpy.isfinite(list(gamma))==True).nonzero()
2799 if len(indxs) >0:
2800 if numpy.nanmean(gamma) > 0.07:
2801 maxp0 = numpy.argmax(signal0*gamma)
2802 maxp1 = numpy.argmax(signal1*gamma)
2803 #print('usa gamma',numpy.nanmean(gamma))
2804 else:
2829 maxp0 = numpy.argmax(signal0)
2805 maxp0 = numpy.argmax(signal0)
2830 maxp1 = numpy.argmax(signal1)
2806 maxp1 = numpy.argmax(signal1)
2831 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
2807 jvelr[h] = (absc[maxp0]+absc[maxp1])/2.
2808 else: jvelr[h] = absc[0]
2832 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
2809 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
2833 #print(maxp0,absc[maxp0],snr0,jvelr[h])
2810 #print(maxp0,absc[maxp0],snr0,jvelr[h])
2834
2811
2835 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
2812 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
2836 else: fd0 = numpy.nan
2813 else: fd0 = numpy.nan
2837
2838 for h in range(nHeights):
2814 for h in range(nHeights):
2839 d = data[:,h]
2815 d = data[:,h]
2840 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
2816 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
@@ -2859,11 +2835,8 class SpectralFitting(Operation):
2859 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
2835 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
2860 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
2836 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
2861 ind += 1
2837 ind += 1
2862 #diagD = numpy.zeros(256)
2838 diagD = numpy.zeros(256)
2863 #if h == 17 :
2839
2864 # for ii in range(256): diagD[ii] = D[ii,ii]
2865 #Dinv=numpy.linalg.inv(D)
2866 #L=numpy.linalg.cholesky(Dinv)
2867 try:
2840 try:
2868 Dinv=numpy.linalg.inv(D)
2841 Dinv=numpy.linalg.inv(D)
2869 L=numpy.linalg.cholesky(Dinv)
2842 L=numpy.linalg.cholesky(Dinv)
@@ -2873,14 +2846,18 class SpectralFitting(Operation):
2873 LT=L.T
2846 LT=L.T
2874
2847
2875 dp = numpy.dot(LT,d)
2848 dp = numpy.dot(LT,d)
2876
2877 #Initial values
2849 #Initial values
2878 data_spc = dataOut.data_spc[coord,:,h]
2850 data_spc = dataOut.data_spc[coord,:,h]
2851 w = data_spc/data_spc
2852 if filec != None:
2853 w = self.weightf.weightfit(w,tini.tm_year,tini.tm_yday,index,h,i)
2879 if (h>6)and(error1[3]<25):
2854 if (h>6)and(error1[3]<25):
2880 p0 = dataOut.data_param[i,:,h-1]
2855 p0 = dataOut.data_param[i,:,h-1]
2881 else:
2856 else:
2882 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))# sin el i(data_spc, constants, i)
2857 p0 = numpy.array(self.library.initialValuesFunction(data_spc*w, constants))# sin el i(data_spc, constants, i)
2883 p0[3] = fd0
2858 p0[3] = fd0
2859 if filec != None:
2860 p0 = self.weightf.Vrfit(p0,tini.tm_year,tini.tm_yday,index,h,i)
2884 try:
2861 try:
2885 #Least Squares
2862 #Least Squares
2886 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
2863 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
@@ -2924,11 +2901,11 class SpectralFitting(Operation):
2924 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
2901 dataOut.data_paramC[0+4*i,ht,0] = numpy.sum((signalpn0/val0_npoints))/n0
2925 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
2902 dataOut.data_paramC[1+4*i,ht,0] = numpy.sum((signalpn1/val1_npoints))/n1
2926
2903
2927 signal0 = (signalpn0-n0) # > 0
2904 signal0 = (signalpn0-n0)
2928 vali = (signal0 < 0).nonzero()
2905 vali = (signal0 < 0).nonzero()
2929 vali = vali[0]
2906 vali = vali[0]
2930 if len(vali) > 0 : signal0[vali] = 0
2907 if len(vali) > 0 : signal0[vali] = 0
2931 signal1 = (signalpn1-n1) #> 0
2908 signal1 = (signalpn1-n1)
2932 vali = (signal1 < 0).nonzero()
2909 vali = (signal1 < 0).nonzero()
2933 vali = vali[0]
2910 vali = vali[0]
2934 if len(vali) > 0 : signal1[vali] = 0
2911 if len(vali) > 0 : signal1[vali] = 0
@@ -2947,6 +2924,7 class SpectralFitting(Operation):
2947
2924
2948 dataOut.data_spc = jspectra
2925 dataOut.data_spc = jspectra
2949 dataOut.spc_noise = my_noises*nProf*M
2926 dataOut.spc_noise = my_noises*nProf*M
2927 if numpy.any(proc): dataOut.spc_noise = my_noises*nProf*M
2950 if getSNR:
2928 if getSNR:
2951 listChannels = groupArray.reshape((groupArray.size))
2929 listChannels = groupArray.reshape((groupArray.size))
2952 listChannels.sort()
2930 listChannels.sort()
@@ -3468,8 +3446,8 class WindProfiler(Operation):
3468
3446
3469 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
3447 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
3470
3448
3471 param = dataOut.data_param
3449 param = dataOut.moments
3472 if dataOut.abscissaList != None:
3450 if numpy.any(dataOut.abscissaList):
3473 absc = dataOut.abscissaList[:-1]
3451 absc = dataOut.abscissaList[:-1]
3474 # noise = dataOut.noise
3452 # noise = dataOut.noise
3475 heightList = dataOut.heightList
3453 heightList = dataOut.heightList
@@ -3610,7 +3588,7 class WindProfiler(Operation):
3610 dataOut.flagNoData = False
3588 dataOut.flagNoData = False
3611 self.__buffer = None
3589 self.__buffer = None
3612
3590
3613 return
3591 return dataOut
3614
3592
3615 class EWDriftsEstimation(Operation):
3593 class EWDriftsEstimation(Operation):
3616
3594
@@ -3638,18 +3616,13 class EWDriftsEstimation(Operation):
3638 y1=y1[vali]
3616 y1=y1[vali]
3639 x = x[vali]
3617 x = x[vali]
3640 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
3618 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
3641
3642 #heiRang1 = x*math.cos(phi[maxid])
3643 x1 = heiRang1
3619 x1 = heiRang1
3644 y11 = f1(x1)
3620 y11 = f1(x1)
3645
3646 y2 = SNR[i,:]
3621 y2 = SNR[i,:]
3647 #print 'snr ', y2
3648 x = heiRang*math.cos(phi[i])
3622 x = heiRang*math.cos(phi[i])
3649 vali= (y2 != -1).nonzero()
3623 vali= (y2 != -1).nonzero()
3650 y2 = y2[vali]
3624 y2 = y2[vali]
3651 x = x[vali]
3625 x = x[vali]
3652 #print 'snr ',y2
3653 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
3626 f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
3654 y21 = f2(x1)
3627 y21 = f2(x1)
3655
3628
@@ -3663,12 +3636,10 class EWDriftsEstimation(Operation):
3663 heiRang = dataOut.heightList
3636 heiRang = dataOut.heightList
3664 velRadial = dataOut.data_param[:,3,:]
3637 velRadial = dataOut.data_param[:,3,:]
3665 velRadialm = dataOut.data_param[:,2:4,:]*-1
3638 velRadialm = dataOut.data_param[:,2:4,:]*-1
3666
3667 rbufc=dataOut.data_paramC[:,:,0]
3639 rbufc=dataOut.data_paramC[:,:,0]
3668 ebufc=dataOut.data_paramC[:,:,1]
3640 ebufc=dataOut.data_paramC[:,:,1]
3669 SNR = dataOut.data_snr
3641 SNR = dataOut.data_snr
3670 velRerr = dataOut.data_error[:,4,:]
3642 velRerr = dataOut.data_error[:,4,:]
3671 #moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
3672 channels = dataOut.channelList
3643 channels = dataOut.channelList
3673 nChan = len(channels)
3644 nChan = len(channels)
3674 my_nbeams = nChan/2
3645 my_nbeams = nChan/2
@@ -3683,7 +3654,6 class EWDriftsEstimation(Operation):
3683 p_w1C = rbufc[1,:]
3654 p_w1C = rbufc[1,:]
3684 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
3655 w_wC = rbufc[2,:]*-1 #*radial_sign(radial EQ 1)
3685 t_wC = rbufc[3,:]
3656 t_wC = rbufc[3,:]
3686 #my_nbeams = 2
3687 if my_nbeams == 1:
3657 if my_nbeams == 1:
3688 w = velRadial[0,:]
3658 w = velRadial[0,:]
3689 winds = velRadial.copy()
3659 winds = velRadial.copy()
@@ -3714,7 +3684,6 class EWDriftsEstimation(Operation):
3714 if len(bad) > 0 :
3684 if len(bad) > 0 :
3715 w_w[bad] = w_wC[bad]
3685 w_w[bad] = w_wC[bad]
3716 w_w_err[bad]= numpy.nan
3686 w_w_err[bad]= numpy.nan
3717 # if my_nbeams == 2:
3718 smooth_eC=ebufc[4,:]
3687 smooth_eC=ebufc[4,:]
3719 p_e0C = rbufc[4,:]
3688 p_e0C = rbufc[4,:]
3720 p_e1C = rbufc[5,:]
3689 p_e1C = rbufc[5,:]
@@ -3745,14 +3714,12 class EWDriftsEstimation(Operation):
3745
3714
3746 hei_aver0 = 218
3715 hei_aver0 = 218
3747 jrange = 450 #900 para HA drifts
3716 jrange = 450 #900 para HA drifts
3748 deltah = 15 # para EWDrifts 15.0 MP 15 #dataOut.spacing(0)
3717 deltah = 15.0 #dataOut.spacing(0) 25 HAD
3749 h0 = 0.0 #dataOut.first_height(0)
3718 h0 = 0.0 #dataOut.first_height(0)
3750 heights = dataOut.heightList
3719 heights = dataOut.heightList
3751 nhei = len(heights)
3720 nhei = len(heights)
3752
3721
3753 range1 = numpy.arange(nhei) * deltah + h0
3722 range1 = numpy.arange(nhei) * deltah + h0
3754
3755 #jhei = WHERE(range1 GE hei_aver0 , jcount)
3756 jhei = (range1 >= hei_aver0).nonzero()
3723 jhei = (range1 >= hei_aver0).nonzero()
3757 if len(jhei[0]) > 0 :
3724 if len(jhei[0]) > 0 :
3758 h0_index = jhei[0][0] # Initial height for getting averages 218km
3725 h0_index = jhei[0][0] # Initial height for getting averages 218km
@@ -3776,7 +3743,6 class EWDriftsEstimation(Operation):
3776
3743
3777 wA = w[h0_index:h0_index+nhei_avg-0]
3744 wA = w[h0_index:h0_index+nhei_avg-0]
3778 wA_err = w_err[h0_index:h0_index+nhei_avg-0]
3745 wA_err = w_err[h0_index:h0_index+nhei_avg-0]
3779
3780 for i in range(5) :
3746 for i in range(5) :
3781 vals = wA[i*h_avgs:(i+1)*h_avgs-0]
3747 vals = wA[i*h_avgs:(i+1)*h_avgs-0]
3782 errs = wA_err[i*h_avgs:(i+1)*h_avgs-0]
3748 errs = wA_err[i*h_avgs:(i+1)*h_avgs-0]
@@ -3796,7 +3762,6 class EWDriftsEstimation(Operation):
3796 wA = wA[6*h_avgs:nhei_avg-0]
3762 wA = wA[6*h_avgs:nhei_avg-0]
3797 wA_err=wA_err[6*h_avgs:nhei_avg-0]
3763 wA_err=wA_err[6*h_avgs:nhei_avg-0]
3798 if my_nbeams == 2 :
3764 if my_nbeams == 2 :
3799
3800 uA = u[h0_index:h0_index+nhei_avg]
3765 uA = u[h0_index:h0_index+nhei_avg]
3801 uA_err=u_err[h0_index:h0_index+nhei_avg]
3766 uA_err=u_err[h0_index:h0_index+nhei_avg]
3802
3767
@@ -3821,12 +3786,9 class EWDriftsEstimation(Operation):
3821 tini=time.localtime(dataOut.utctime)
3786 tini=time.localtime(dataOut.utctime)
3822 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
3787 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
3823 nfile = fileDrifts+'/jro'+datefile+'drifts.txt'
3788 nfile = fileDrifts+'/jro'+datefile+'drifts.txt'
3824
3825 f1 = open(nfile,'a')
3789 f1 = open(nfile,'a')
3826
3827 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
3790 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
3828 driftavgstr=str(dataOut.drifts_avg)
3791 driftavgstr=str(dataOut.drifts_avg)
3829
3830 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
3792 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
3831 numpy.savetxt(f1,dataOut.drifts_avg,fmt='%10.2f')
3793 numpy.savetxt(f1,dataOut.drifts_avg,fmt='%10.2f')
3832 f1.close()
3794 f1.close()
@@ -5243,3 +5205,137 class IGRFModel(Operation):
5243 mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
5205 mkfact_short_2020.mkfact(dataOut.year,dataOut.h,dataOut.bfm,dataOut.thb,dataOut.bki,dataOut.MAXNRANGENDT)
5244
5206
5245 return dataOut
5207 return dataOut
5208
5209 class MergeProc(ProcessingUnit):
5210
5211 def __init__(self):
5212 ProcessingUnit.__init__(self)
5213
5214 def run(self, attr_data, attr_data_2 = None, attr_data_3 = None, attr_data_4 = None, attr_data_5 = None, mode=0):
5215
5216 self.dataOut = getattr(self, self.inputs[0])
5217 data_inputs = [getattr(self, attr) for attr in self.inputs]
5218 #print(self.inputs)
5219 #print(numpy.shape([getattr(data, attr_data) for data in data_inputs][1]))
5220 #exit(1)
5221 if mode==0:
5222 data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs])
5223 setattr(self.dataOut, attr_data, data)
5224
5225 if mode==1: #Hybrid
5226 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
5227 #setattr(self.dataOut, attr_data, data)
5228 setattr(self.dataOut, 'dataLag_spc', [getattr(data, attr_data) for data in data_inputs][0])
5229 setattr(self.dataOut, 'dataLag_spc_LP', [getattr(data, attr_data) for data in data_inputs][1])
5230 setattr(self.dataOut, 'dataLag_cspc', [getattr(data, attr_data_2) for data in data_inputs][0])
5231 setattr(self.dataOut, 'dataLag_cspc_LP', [getattr(data, attr_data_2) for data in data_inputs][1])
5232 #setattr(self.dataOut, 'nIncohInt', [getattr(data, attr_data_3) for data in data_inputs][0])
5233 #setattr(self.dataOut, 'nIncohInt_LP', [getattr(data, attr_data_3) for data in data_inputs][1])
5234 '''
5235 print(self.dataOut.dataLag_spc_LP.shape)
5236 print(self.dataOut.dataLag_cspc_LP.shape)
5237 exit(1)
5238 '''
5239
5240 #self.dataOut.dataLag_spc_LP = numpy.transpose(self.dataOut.dataLag_spc_LP[0],(2,0,1))
5241 #self.dataOut.dataLag_cspc_LP = numpy.transpose(self.dataOut.dataLag_cspc_LP,(3,1,2,0))
5242 '''
5243 print("Merge")
5244 print(numpy.shape(self.dataOut.dataLag_spc))
5245 print(numpy.shape(self.dataOut.dataLag_spc_LP))
5246 print(numpy.shape(self.dataOut.dataLag_cspc))
5247 print(numpy.shape(self.dataOut.dataLag_cspc_LP))
5248 exit(1)
5249 '''
5250 #print(numpy.sum(self.dataOut.dataLag_spc_LP[2,:,164])/128)
5251 #print(numpy.sum(self.dataOut.dataLag_cspc_LP[0,:,30,1])/128)
5252 #exit(1)
5253 #print(self.dataOut.NDP)
5254 #print(self.dataOut.nNoiseProfiles)
5255
5256 #self.dataOut.nIncohInt_LP = 128
5257 self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
5258 self.dataOut.nIncohInt_LP = self.dataOut.nIncohInt
5259 self.dataOut.NLAG = 16
5260 self.dataOut.NRANGE = 200
5261 self.dataOut.NSCAN = 128
5262 #print(numpy.shape(self.dataOut.data_spc))
5263
5264 #exit(1)
5265
5266 if mode==2: #HAE 2022
5267 data = numpy.sum([getattr(data, attr_data) for data in data_inputs],axis=0)
5268 setattr(self.dataOut, attr_data, data)
5269
5270 self.dataOut.nIncohInt *= 2
5271 #meta = self.dataOut.getFreqRange(1)/1000.
5272 self.dataOut.freqRange = self.dataOut.getFreqRange(1)/1000.
5273
5274 #exit(1)
5275
5276 if mode==4: #Hybrid LP-SSheightProfiles
5277 #data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs],axis=1)
5278 #setattr(self.dataOut, attr_data, data)
5279 setattr(self.dataOut, 'dataLag_spc', getattr(data_inputs[0], attr_data)) #DP
5280 setattr(self.dataOut, 'dataLag_cspc', getattr(data_inputs[0], attr_data_2)) #DP
5281 setattr(self.dataOut, 'dataLag_spc_LP', getattr(data_inputs[1], attr_data_3)) #LP
5282 #setattr(self.dataOut, 'dataLag_cspc_LP', getattr(data_inputs[1], attr_data_4)) #LP
5283 #setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5284 setattr(self.dataOut, 'data_acf', getattr(data_inputs[1], attr_data_5)) #LP
5285 #print("Merge data_acf: ",self.dataOut.data_acf.shape)
5286 #exit(1)
5287 #print(self.dataOut.data_spc_LP.shape)
5288 #print("Exit")
5289 #exit(1)
5290 #setattr(self.dataOut, 'dataLag_cspc', [getattr(data, attr_data_2) for data in data_inputs][0])
5291 #setattr(self.dataOut, 'dataLag_cspc_LP', [getattr(data, attr_data_2) for data in data_inputs][1])
5292 #setattr(self.dataOut, 'nIncohInt', [getattr(data, attr_data_3) for data in data_inputs][0])
5293 #setattr(self.dataOut, 'nIncohInt_LP', [getattr(data, attr_data_3) for data in data_inputs][1])
5294 '''
5295 print(self.dataOut.dataLag_spc_LP.shape)
5296 print(self.dataOut.dataLag_cspc_LP.shape)
5297 exit(1)
5298 '''
5299 '''
5300 print(self.dataOut.dataLag_spc_LP[0,:,100])
5301 print(self.dataOut.dataLag_spc_LP[1,:,100])
5302 exit(1)
5303 '''
5304 #self.dataOut.dataLag_spc_LP = numpy.transpose(self.dataOut.dataLag_spc_LP[0],(2,0,1))
5305 #self.dataOut.dataLag_cspc_LP = numpy.transpose(self.dataOut.dataLag_cspc_LP,(3,1,2,0))
5306 '''
5307 print("Merge")
5308 print(numpy.shape(self.dataOut.dataLag_spc))
5309 print(numpy.shape(self.dataOut.dataLag_spc_LP))
5310 print(numpy.shape(self.dataOut.dataLag_cspc))
5311 print(numpy.shape(self.dataOut.dataLag_cspc_LP))
5312 exit(1)
5313 '''
5314 #print(numpy.sum(self.dataOut.dataLag_spc_LP[2,:,164])/128)
5315 #print(numpy.sum(self.dataOut.dataLag_cspc_LP[0,:,30,1])/128)
5316 #exit(1)
5317 #print(self.dataOut.NDP)
5318 #print(self.dataOut.nNoiseProfiles)
5319
5320 #self.dataOut.nIncohInt_LP = 128
5321 #self.dataOut.nProfiles_LP = 128#self.dataOut.nIncohInt_LP
5322 self.dataOut.nProfiles_LP = 16#28#self.dataOut.nIncohInt_LP
5323 self.dataOut.nProfiles_LP = self.dataOut.data_acf.shape[1]#28#self.dataOut.nIncohInt_LP
5324 self.dataOut.NSCAN = 128
5325 self.dataOut.nIncohInt_LP = self.dataOut.nIncohInt*self.dataOut.NSCAN
5326 #print("sahpi",self.dataOut.nIncohInt_LP)
5327 #exit(1)
5328 self.dataOut.NLAG = 16
5329 self.dataOut.NRANGE = self.dataOut.data_acf.shape[-1]
5330
5331 #print(numpy.shape(self.dataOut.data_spc))
5332
5333 #exit(1)
5334 if mode==5:
5335 data = numpy.concatenate([getattr(data, attr_data) for data in data_inputs])
5336 setattr(self.dataOut, attr_data, data)
5337 data = numpy.concatenate([getattr(data, attr_data_2) for data in data_inputs])
5338 setattr(self.dataOut, attr_data_2, data)
5339 #data = numpy.concatenate([getattr(data, attr_data_3) for data in data_inputs])
5340 #setattr(self.dataOut, attr_data_3, data)
5341 #print(self.dataOut.moments.shape,self.dataOut.data_snr.shape,self.dataOut.heightList.shape) No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now