<jupyter_start><jupyter_text>Routine to get parameters<jupyter_code>TRES = 5#10
VRES = 3#6

def CalSys(year, month, day, hour, TRES=TRES, VRES=VRES, rlim=4, doshow=True, rlimmax=7):

    path_in = "I:/PHD/Lidar/Processing_V3/Signals/"+str(TRES)+"min"+str(VRES)+"bins/"
    filename1 = path_in+"Par90_bc/Par90bc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
    filename2 = path_in+"Par10_bc/Par10bc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
    filename3 = path_in+"Per_bc/Perbc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"

    filenameSNR = path_in+"SNR/SNR_Par90_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
    
    r   = np.loadtxt(path_in+"R_"+str(TRES)+"min_"+str(VRES)+"bins.dat")
    Mat_Spar = np.loadtxt(filename1)
    Mat_Spar10 = np.loadtxt(filename2)
    Mat_Sper = np.loadtxt(filename3)
    SNR_Spar = np.loadtxt(filenameSNR)

    pix_r = np.squeeze(np.where((r >= rlim)&(r <= rlimmax)))[0]
    
    #SNR_par2 = np.ma.masked_where(~(SNR_par>30),SNR_par)
    #pix_r2 = np.squeeze(np.where(r >= 0.3))[0]
    
    Spar = Mat_Spar[int((60./TRES)*hour),:]
    Spar10 = Mat_Spar10[int((60./TRES)*hour),:]
    Sper = Mat_Sper[int((60./TRES)*hour),:]
    
    SNR_Spar2 = SNR_Spar[int((60./TRES)*hour),:]
    
    #pixK = np.where(~(SNR_Spar2>30))
    
    r_sm = BG_corr(Spar,r, year, month, day,rcf0 = 8)[2]
    sm = BG_corr(Spar,r, year, month, day,rcf0 = 8)[3]

    f_sm = interp1d(r_sm, sm, fill_value='extrapolate')

    sm_lidar = f_sm(r[int(11*3./VRES):])#/100000
    sm_lidar[np.isnan(sm_lidar)] = 0
    bm = sm_lidar/((8*np.pi/3)*1.0401)
    bm_par = 0.996*bm
    bm_per =0.00366*bm_par
    T2 = np.exp((-2*np.cumsum(sm_lidar)*(r[1]-r[0])))  
    

    Cpar_ = Spar[int(11*3./VRES):]/(bm_par*T2)
    Cper_ = Sper[int(11*3./VRES):]/(bm_per*T2)
    
    pix = np.squeeze(np.where((r[int(11*3./VRES):]>rlim) & (r[int(11*3./VRES):]<rlimmax) & (SNR_Spar2[int(11*3./VRES):]>10)))
    
    #print pix, Cpar_
    
    Cpar = np.nanmedian(Cpar_[pix])#np.nanmean(Cpar_[pix_r:])
    Cper = np.nanmedian(Cper_[pix])# np.nanmean(Cper_[pix_r:])
    Cpar_std = np.nanstd(Cpar_[pix])#np.nanstd(Cpar_[pix_r:])
    Cper_std = np.nanstd(Cper_[pix])#np.nanstd(Cper_[pix_r:])  
    
    K = np.nanmedian((0.00366/(Sper/Spar)))#[pixK])
    #pylab.plot(r[pixK],(0.00366/(Sper/Spar))[pixK])
    #pylab.show()

    pix = np.where((r>0.3) & (SNR_Spar2>30))
    CC = (0.00366/(Sper/Spar))[pix]  
    
    K2 = np.nanmedian(CC)
    
    if doshow == True:
        
        pylab.plot(r[int(11*3./VRES):],1/(bm_par*T2))
        pylab.plot(r[int(11*3./VRES):],1/(bm_per*T2))
        #pylab.axis([3,10,1,0.01])
        pylab.show()
        
        pylab.plot(r[:],Spar)
        pylab.plot(r[:],Sper)
        pylab.axis([3,10,1,1000])
        pylab.xlabel("Range [km]")
        pylab.ylabel("r2-corrected Signal")
        pylab.show()

        pylab.plot(r[:],0.00366/(((Sper))/(Spar)))
        pylab.axis([3,10,0,1])
        pylab.xlabel("Range [km]")
        pylab.ylabel("Dep.Ratio")
        pylab.show()

        pylab.plot(r[int(11*3./VRES):][pix_r:],Cpar_[pix_r:])
        pylab.xlabel("Range [km]")
        pylab.ylabel("Cpar")
        pylab.show()

        pylab.plot(r[int(11*3./VRES):][pix_r:],Cper_[pix_r:])
        pylab.xlabel("Range [km]")
        pylab.ylabel("Cper")
        pylab.show()    
    
        print "Cpar = ", Cpar,'/std = ',Cpar_std,"Cper = ",Cper,'/std = ',Cper_std,"Cpar/Cper = ", Cpar/Cper
    return [[Cpar,Cper], [Spar,Spar10, Sper], r, K, SNR_Spar2, K2, [r[int(11*3./VRES):],sm_lidar]]


dates = np.zeros(shape = (29,4)) 

dates[0,:] = np.array([2017,2,25,14]) #[YYYY,M,D,H]

dates[1,:] = np.array([2017,3,8,14]) #[YYYY,M,D,H]
dates[2,:] = np.array([2017,3,18,13]) #[YYYY,M,D,H]
dates[3,:] = np.array([2017,3,20,16]) #[YYYY,M,D,H]
dates[4,:] = np.array([2017,3,21,14]) #[YYYY,M,D,H]
dates[5,:] = np.array([2017,3,24,12]) #[YYYY,M,D,H]

dates[6,:] = np.array([2017,4,6,10]) #[YYYY,M,D,H]
dates[7,:] = np.array([2017,4,9,19]) #[YYYY,M,D,H]
#dates[7,:] = np.array([2017,4,10,12]) #[YYYY,M,D,H]
dates[8,:] = np.array([2017,4,11,11]) #[YYYY,M,D,H]
dates[9,:] = np.array([2017,4,18,13]) #[YYYY,M,D,H]

dates[10,:] = np.array([2017,5,2,8]) #[YYYY,M,D,H]
dates[11,:] = np.array([2017,5,3,8]) #[YYYY,M,D,H]
dates[12,:] = np.array([2017,5,4,10]) #[YYYY,M,D,H]
dates[13,:] = np.array([2017,5,7,12]) #[YYYY,M,D,H]
dates[14,:] = np.array([2017,5,11,13]) #[YYYY,M,D,H]
dates[15,:] = np.array([2017,5,22,18]) #[YYYY,M,D,H]
dates[16,:] = np.array([2017,5,23,8]) #[YYYY,M,D,H]

dates[17,:] = np.array([2017,6,1,14]) #[YYYY,M,D,H]
dates[18,:] = np.array([2017,6,6,8]) #[YYYY,M,D,H]
dates[19,:] = np.array([2017,6,11,11]) #[YYYY,M,D,H]
dates[20,:] = np.array([2017,6,12,18]) #[YYYY,M,D,H]
dates[21,:] = np.array([2017,6,13,14]) #[YYYY,M,D,H]
dates[22,:] = np.array([2017,6,18,18]) #[YYYY,M,D,H]

dates[23,:] = np.array([2017,7,9,14]) #[YYYY,M,D,H]
dates[24,:] = np.array([2017,7,10,16]) #[YYYY,M,D,H]
dates[25,:] = np.array([2017,7,18,15]) #[YYYY,M,D,H]

dates[26,:] = np.array([2017,8,9,10]) #[YYYY,M,D,H]
dates[27,:] = np.array([2017,8,14,14]) #[YYYY,M,D,H]
dates[28,:] = np.array([2017,8,22,10]) #[YYYY,M,D,H]

Cs = np.zeros(shape = (np.shape(dates)[0], 2))
CCs = np.zeros(shape = (np.shape(dates)[0]))
Ks = np.zeros(shape = (np.shape(dates)[0]))

epoch_time = []

for i in range(np.shape(dates)[0]):
    Cs[i,:],P,r,Ks[i], SNR, CCs[i], sigma = CalSys(int(dates[i,0]), int(dates[i,1]), int(dates[i,2]),int(dates[i,3]), TRES = TRES, VRES = VRES, doshow = False)
    utc_time = time.strptime(str(int(dates[i,0]))+"-"+str(int(dates[i,1])).zfill(2)+"-"+str(int(dates[i,2])).zfill(2)+" "+str(int(dates[i,3])).zfill(2)+":00:00", "%Y-%m-%d %H:%M:%S")
    epoch_time.append(timegm(utc_time))
epoch_time = np.array(epoch_time)

doy = 31+dates[0,2]+dates[0,3]/24.+(epoch_time-epoch_time[0])/(3600.*24)<jupyter_output>C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:57: RuntimeWarning: divide by zero encountered in divide
C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:62: RuntimeWarning: divide by zero encountered in divide
C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:48: RuntimeWarning: invalid value encountered in greater
C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:57: RuntimeWarning: invalid value encountered in divide
C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:61: RuntimeWarning: invalid value encountered in greater
C:\Users\duran\AppData\Local\Continuum\Anaconda2\lib\site-packages\ipykernel\__main__.py:62: RuntimeWarning: invalid value encountered in divide<jupyter_text>Save Parameters<jupyter_code>filename = "I:/PHD/Lidar/Processing_V3/Calibration_sys/CalSys_"+str(TRES)+"min_"+str(VRES)+"bins_1.csv"
if os.path.isfile(filename) == False: 
    txt = open(filename,"w")
    txt.write('timestamp,EPOCH_time,DOY,Cpar,Cper,K'+"\n")
    for i in range(np.shape(dates)[0]):
        txt.write(str(int(dates[i,0]))+"-"+str(int(dates[i,1])).zfill(2)+"-"+str(int(dates[i,2])).zfill(2)+" "+str(int(dates[i,3])).zfill(2)+":00:00"+","+str(epoch_time[i])+","+str(doy[i])+","+str(Cs[i,0])+","+str(Cs[i,1])+","+str(CCs[i])+"\n")
    txt.close()    
# PLots
fig, ax1 = pylab.subplots(figsize=(13,8))
ax1.plot(doy,Cs[:,0], "-o", label = r'$C_{\parallel}$')
ax1.plot(doy,Cs[:,1], "-o", label = r'$C_{\perp}$')
ax1.set_xlabel("Day of the year 2017")
ax1.set_ylabel(r'$C_{\parallel,\perp}$')
#pylab.show()
pylab.axis([50,250,0,4.5e10])
ax2 = ax1.twinx()
ax2.plot(doy,Cs[:,0], "-o", label = "C"+r'$_{\parallel}$', color = "blue")
ax2.plot(doy,Cs[:,1], "-o", label = "C"+r'$_{\perp}$', color = "green")
ax2.plot(doy, CCs, "-", color = 'black', label = r'$K$')
pylab.axis([50,250,0,0.5])
ax2.set_ylabel(r'$K = C_{\parallel}/C_{\perp}$')
ax2.legend(frameon = False)
pylab.title("System Calibration")
pylab.savefig("plots/System_calibration"+str(TRES)+"min_"+str(VRES)+"bins_3.png",format="png",bbox_inches = 'tight', dpi=900)
pylab.show()<jupyter_output><empty_output><jupyter_text>Following codes are only for tests<jupyter_code>#Overlapping correction
##### output Temporal RESolution
TRES = 10 #min. 

##### output Temporal resolution
VRES = 6 #bins, 1bin = 3.8m # try only 1,2,3,6 bins

year = 2017
month = 3
day  = 8
hour = 14
path_in = "I:/PHD/Lidar/Processing_V3/Signals/"+str(TRES)+"min"+str(VRES)+"bins/"

filename = path_in+"Par90_bc/Par90bc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
Par = np.loadtxt(filename)

Spar = Par[60/TRES*14,:]

r   = np.loadtxt(path_in+"R_"+str(TRES)+"min_"+str(VRES)+"bins.dat")

r_sm = BG_corr(Spar,r, year, month, day,rcf0 = 8)[2]
sm = BG_corr(Spar,r, year, month, day,rcf0 = 8)[3]

f_sm = interp1d(r_sm, sm, fill_value='extrapolate')

sm_lidar = f_sm(r[int(11*3./VRES):])#/100000

C,P,r,K, SNR, CC, sigma = CalSys(year, month, day,hour, TRES = TRES, VRES = VRES, doshow = False)

B = P[0][int(11*3./VRES):]/C[0]
Bm = 0.996*sm_lidar/((8*np.pi/3)*1.0401)*np.exp(-2*np.cumsum(sm_lidar)*(r[1]-r[0]))

#O = np.loadtxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"_08032017_14UTC.dat")
O3 = np.loadtxt("C:/Users/duran/Documents/PhD/CODES/Python27/Notebooks/LIDAR/OverlapCorrBr_15min-6pix.txt")
r3 = np.loadtxt("C:/Users/duran/Documents/PhD/CODES/Python27/Notebooks/LIDAR/r_15min-6pix.txt")

f_0 = interp1d(r3[3:]/1000.,O3, fill_value='extrapolate')
O4 =f_0(r[int(11*3./VRES):])

pylab.plot(r[int(11*3./VRES):],Bm)
pylab.plot(r[int(11*3./VRES):],B)
line_Bm = np.exp(lm[0]*r[int(11*3./VRES):])*np.exp(lm[1]) 
pylab.plot(r[int(11*3./VRES):],line_Bm)
pylab.yscale("log")
pylab.show()

O5 = (line_Bm/B)
O2 = (Bm/B)

#pix_lin = np.where((r[int(11*3./VRES):]>4) & (r[int(11*3./VRES):]<7))

#X = r[int(11*3./VRES):][pix_lin]
#Y = np.log(B[pix_lin])

x = r[int(11*3./VRES):]
p = np.polyfit(1/x**2, O5, 2)
O5_fit = p[0]/x**4+p[1]/x**2+p[2]

p = np.polyfit(1/x**2, O2, 2)
O2_fit = p[0]/x**4+p[1]/x**2+p[2]

x = r[int(11*3./VRES):]
p = np.polyfit(1/x**2, O5, 2)
O5_fit = p[0]/x**4+p[1]/x**2+p[2]

p = np.polyfit(1/x**2, O4, 2)
O4_fit = p[0]/x**4+p[1]/x**2+p[2]

#pylab.plot(1/x**2,O2)
#pylab.plot(1/x**2,O2_fit)
#pylab.plot(x,O2)
#pylab.plot(x,O2_fit)
#pylab.show()

#pylab.plot(1/x**2,O5)
#pylab.plot(1/x**2,O5_fit)
#pylab.show()

#pylab.plot(r[int(11*3./VRES):],O2)
#pylab.plot(r[int(11*3./VRES):],O2_fit)
pylab.plot(r3[3:]/1000.,O3)
pylab.plot(x,O4_fit)
pylab.plot(x,O2_fit)
#pylab.plot(x,O5)
pylab.plot(x,O5_fit)
pylab.axis([0,2,0,50])
pylab.show()

np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC.dat",O2_fit)
np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC_r.dat",r[int(11*3./VRES):])

np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC2.dat",O4_fit)
np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC_r2.dat",r[int(11*3./VRES):])

np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC3.dat",O5_fit)
np.savetxt(path_in+"Overlapping_"+str(TRES)+"min_"+str(VRES)+"bins_08032017_14UTC_r3.dat",r[int(11*3./VRES):])<jupyter_output><empty_output><jupyter_text>Cross-talk correction<jupyter_code>Par =  np.loadtxt("I:/PHD/Lidar/Processing_V3/SysCalibrated/10min6bins/Bpar/Bpar_sc_O_20170226_10min_6bins.dat")
Per =  np.loadtxt("I:/PHD/Lidar/Processing_V3/SysCalibrated/10min6bins/Bper/Bper_sc_O_20170226_10min_6bins.dat")
SNR =  np.loadtxt("I:/PHD/Lidar/Processing_V3/Signals/10min6bins/SNR/SNR_Par10_20170226_10min_6bins.dat")[:,5:]
r = np.loadtxt("I:/PHD/Lidar/Processing_V3/Signals/10min6bins/R_10min_6bins.dat")[5:]
dr = Per/Par
pix
SNR[3*8,20:39]
np.squeeze(np.where(r <= pix[-1]))[-1]
h = 8

pix  = np.squeeze(np.where(SNR[6*h,:]> 25))
    
r0 = np.squeeze(np.where(r >= r[pix[0]]))[0]
r1 = np.squeeze(np.where(r <= r[pix[-1]]))[-1]

fig, ax1 = pylab.subplots(figsize=(13,8))    
pylab.plot(r[r0:r1],1000*Par[6*h,r0:r1])
pylab.plot(r[r0:r1],1000*40*Per[6*h,r0:r1])
pylab.plot(r[r0:r1],dr[6*h,r0:r1])
pylab.axis([r[pix[0]],r[pix[-1]],0,0.12])
pylab.show()

fig, ax1 = pylab.subplots(figsize=(13,8))
pylab.plot(Par[6*h,r0:r1],40*Per[6*h,r0:r1], "o")

pylab.show()


pylab.plot(r,SNR)
pylab.plot(r,SNR*0+30)

pylab.axis([1,10,0,300])
pylab.show()

np.where((r[int(10*3./VRES):]>5) & (SNR[int(10*3./VRES):]>30))
path_in = "I:/PHD/Lidar/Processing_V3/Signals/"+str(TRES)+"min"+str(VRES)+"bins/"
filename1 = path_in+"Par90_bc/Par90bc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
filename3 = path_in+"Per_bc/Perbc_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"
filenameSNR = path_in+"SNR/SNR_Par90_"+str(year)+str(month).zfill(2)+str(day).zfill(2)+"_"+str(TRES)+"min_"+str(VRES)+"bins.dat"

r   = np.loadtxt(path_in+"R_"+str(TRES)+"min_"+str(VRES)+"bins.dat")
Mat_Spar = np.loadtxt(filename1)
Mat_Sper = np.loadtxt(filename3)
SNR_par = np.loadtxt(filenameSNR)
hours = np.linspace(0,24,np.shape(Mat_Spar)[0])

pylab.pcolor(hours,r,K*np.transpose(Mat_Sper/Mat_Spar), vmin = 0, vmax = 0.1)
pylab.show()
pylab.pcolor(hours,r,np.transpose(SNR_par),vmin = 0, vmax = 10000)
pylab.show()
SNR_par2 = np.ma.masked_where(~((SNR_par>30) & (r>0.3)),SNR_par)
#SNR_par2 = np.ma.masked_where(~(((SNR_par2) < 0) & ((SNR_par2) >= 0)),SNR_par2)

pylab.pcolor(hours,r,np.transpose(SNR_par2),vmin = 0, vmax = 100)
pylab.colorbar()
pylab.show()


print Cal1[0]
print Cal2[0]
print Cal3[0]
print Cal1[0][0]/Cal1[0][1],Cal2[0][0]/Cal2[0][1],Cal3[0][0]/Cal3[0][1]

pylab.plot(Cal1[2][2],Cal1[2][0]/Cal1[0][0])
pylab.plot(Cal2[2][2],Cal2[2][0]/Cal2[0][0])
pylab.plot(Cal3[2][2],Cal3[2][0]/Cal3[0][0])
pylab.show()

pylab.plot(Cal1[2][2],Cal1[2][1]/Cal1[0][1])
pylab.plot(Cal2[2][2],Cal2[2][1]/Cal2[0][1])
pylab.plot(Cal3[2][2],Cal3[2][1]/Cal3[0][1])
pylab.show()

pylab.plot(Cal1[2][2],(Cal1[2][1]/Cal1[2][0])*(Cal1[0][0]/Cal1[0][1]))
pylab.plot(Cal2[2][2],(Cal2[2][1]/Cal2[2][0])*(Cal2[0][0]/Cal2[0][1]))
pylab.plot(Cal3[2][2],(Cal3[2][1]/Cal3[2][0])*(Cal3[0][0]/Cal3[0][1]))
pylab.axis([1,10,0,0.01])
pylab.show()

K1 = np.nanmedian((0.00366/(Cal1[2][1]/Cal1[2][0])))
K2 = np.nanmedian((0.00366/(Cal2[2][1]/Cal2[2][0])))
K3 = np.nanmedian((0.00366/(Cal3[2][1]/Cal3[2][0])))


pylab.plot(Cal1[2][2],K1*(Cal1[2][1]/Cal1[2][0]))
pylab.plot(Cal2[2][2],K2*(Cal2[2][1]/Cal2[2][0]))
pylab.plot(Cal3[2][2],K3*(Cal3[2][1]/Cal3[2][0]))
pylab.axis([1,10,0,0.05])
pylab.show()
Cpar =  331077.319291 /std =  21207.0732497 Cper =  2428318.65784 /std =  816402.273351 Cpar/Cper =  0.136340145566
Cpar =  337426.257255 /std =  9169.08627914 Cper =  2197662.30654 /std =  373970.868415 Cpar/Cper =  0.153538719871<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.25976549413735345, 'fraction_numerical': 0.09460636515912897, 'mean_word_length': 5.884686346863469, 'pattern_counts': {'":': 2, '<': 22, '<?xml version=': 0, '>': 32, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '330308', 'license': 'MIT', 'n_tokens_mistral': 8163, 'n_tokens_neox': 7566, 'n_words': 699, 'path': 'Notebooks/System_calibration.ipynb', 'repo_name': 'cduranalarcon/Lidar_DDU', 'size': '14925'}