Banda de 60m Agosto 2019
import datetime as dt
import matplotlib.pyplot as plt
import sys
import numpy as np
station=sys.argv[1]
print(station)
WSPR_file=open("/home/andres/WSPR/60m_dp0gvn/"+station,'r')
date_list=[]
hour_list=[]
SNR_list=[]
for l in WSPR_file:
line=l.split()
date_list.append(dt.datetime.strptime(line[0]+' '+line[1],"%Y-%m-%d %H:%M"))#convert string date to datetime format
hour_list.append(line[1])
SNR_list.append(line[4])
init_date=min(date_list) #search min and max dates of the list of reports, from datetime list of objects
print(init_date)
end_date=max(date_list)
print(end_date)
h_list=[]
hourly=[]
init_time=dt.datetime.strptime("00:00","%H:%M")
for i in hour_list:
time=dt.datetime.strptime(i,"%H:%M")
h_list.append(((time-init_time).total_seconds())/3600)
h=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]
h_m=[0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5,22.5,23.5]
h_m_ext=[0,0.5,1.5,2.5,3.5,4.5,5.5,6.5,7.5,8.5,9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5,17.5,18.5,19.5,20.5,21.5,22.5,23.5,24]
mean=[]
mean_ext=[]
stdev=[]
stdev_ext=[]
n_decodes=[]
min_signal=[]
max_signal=[]
min_signal_ext=[]
max_signal_ext=[]
color_decodes=[]
freq=[]
for i in range(0,24):
WSPR_file2=open("/home/andres/WSPR/60m_dp0gvn/"+station,'r')
signal=[]
for l in WSPR_file2:
a=l.split()
freq.append(float(a[3]))
if(dt.datetime.strptime(a[0]+' ' +a[1],"%Y-%m-%d %H:%M").hour==i):
signal.append(int(a[4]))
if (len(signal)!=0):
mu=np.mean(signal)
dev=np.std(signal)
min_sig=mu-dev
max_sig=mu+dev
n=len(signal)
print(str(i)+":"+str(len(signal)))
else:
mu=-30.0
dev=0.0
min_sig=mu-dev
max_sig=mu+dev
print(str(i)+":"+"none")
n=0
mean.append(mu)
stdev.append(dev)
min_signal.append(min_sig)
max_signal.append(max_sig)
n_decodes.append(n)
WSPR_file.close()
m=(mean[0]+mean[-1])*0.5
mean_ext.append(m)
for i in range(0,len(mean)):
mean_ext.append(mean[i])
mean_ext.append(m)
m=((mean[0]+mean[-1])*0.5)-((stdev[0]+stdev[-1])*0.5)
min_signal_ext.append(m)
for i in range(0,len(min_signal)):
min_signal_ext.append(min_signal[i])
min_signal_ext.append(m)
m=((mean[0]+mean[-1])*0.5)+((stdev[0]+stdev[-1])*0.5)
max_signal_ext.append(m)
for i in range(0,len(max_signal)):
max_signal_ext.append(max_signal[i])
max_signal_ext.append(m)
mean_freq=np.mean(freq)
max_decodes=max(n_decodes)
#Color represents porcentage of spots from total
#If lot of spots green, else degrades through yellow, orange and red
for i in n_decodes:
if (float(i)>=max_decodes*0.75):
color_decodes.append("g")
elif (max_decodes*0.5<=float(i) and float(i)<max_decodes*0.75):
color_decodes.append("y")
#elif (max_decodes*0.25<=float(i) and float(i)<max_decodes*0.5):
#color_decodes.append("orange")
else:
color_decodes.append("r")
print(color_decodes)
x=h_list
y=SNR_list
y_int=[]
x_min=0
x_max=24
for i in y: #generate list in int format for SNR levels
y_int.append(int(i))
y_min=min(y_int) #Only works min() if integer
y_max=max(y_int)
fig, ax = plt.subplots(ncols=1, sharey=True, figsize=(7, 4))
fig.subplots_adjust(hspace=0.5, left=0.07, right=0.97, top=0.9)
#ax.errorbar(h_m,mean,xerr=0.5,yerr=stdev,fmt="o--",ecolor=color_decodes)
ax.errorbar(h_m_ext,mean_ext,xerr=0,yerr=0,fmt="-.")
ax.grid(linestyle='--',alpha=0.5)
#ax.fill_between(h_m,min_signal,max_signal,alpha=0.15,color='k')
#????????????????
#min_signal.append(-25)
#max_signal.append(-15)
#escalonado
#ax.fill_between(h_m_ext,min_signal_ext,max_signal_ext,alpha=0.15,color='k',step="mid")
#Interpolado linealmente
ax.fill_between(h_m_ext,min_signal_ext,max_signal_ext,alpha=0.15,color='k')
ax.scatter(h_m,mean,color=color_decodes,marker='o')
#ax.plot(min_signal)
#ax.plot(max_signal)
#Axis texts
ax.set_xlim(left=0,right=24)
ax.set_xlabel("UTC time\nVertical lines are local sunrise/sunset, Red:"+station+", Black:LU3HO")
ax.set_ylabel("Signal to Noise Ratio (dB)")
ax.set_title("WSPR Analysis of LU3HO received by "+station+" from "+str(init_date)+" to "+str(end_date)+"\nShadow is +/- sigma(SNR) - Dotted blue line is mean(SNR) - Dot is mean, if green: enough data, yellow: little data, red:insufficient data\nMean=-30dB is plot when there is no data - Numbers above mean dots are number of spots")
#Print above each mean point the number of received/transmitted spots.
for i in h:
print(i,n_decodes[i])
plt.text(i+0.35,mean[i]+0.5,n_decodes[i])
#Annotate Mean WSPR frequency on top of graph
ax.annotate("%.4f MHz"%mean_freq,
xy=(0, -30), xycoords='data',
xytext=(0.9, 0.95), textcoords='axes fraction',
horizontalalignment='right', verticalalignment='top',fontsize=17)
#Hour axis
ax.set_xticks([0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24])
#SNR Axis
ax.set_yticks([-30,-25,-20,-15])
#Sunrise/Sunset vertical lines....put a better approach
ax.axvline(9,color='r',linestyle='-.')
ax.axvline(16.25,color='r',linestyle='-.')
ax.axvline(11.2,color='k',linestyle='-.')
ax.axvline(22.2,color='k',linestyle='-.')
plt.show()