Skip to content
Snippets Groups Projects
Commit e13112d0 authored by Carlo Guidi's avatar Carlo Guidi
Browse files

Add acoustics monitoring

parent 1bc445dc
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 271 additions and 0 deletions
scripts/acoustics.py 0 → 100644
+ 271
− 0
View file @ e13112d0
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import colors
import km3pipe as kp
import time
import os
from datetime import datetime
import argparse
parser = argparse.ArgumentParser(description='Online_monitoring')
parser.add_argument('-d', dest='det_id', type=str, nargs=1, required=True,
help='The detector ID (e.g. D0ARCA001)')
parser.add_argument('-o', dest='plot_dir', type=str, nargs=1, required=True,
help='Directory in which the plot is saved')
args = parser.parse_args()
detid = args.det_id[0]
directory = args.plot_dir[0]
def diff(first, second):
second = set(second)
return [item for item in first if item not in second]
#detid='D_ORCA005'
#print(time.time())
db = kp.db.DBManager()
sds = kp.db.StreamDS()
table=db.run_table(detid)
AB=[12,14,16] # Acoustic Beacons
DOM=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18] # DOMs
DU=[1,2,3,4,5] # DUs
#AB=[16] # Acoustic Beacons
#DOM=[7] #
#DU=[3] # DUs
check=True
while check==True:
minrun=table["RUN"][len(table["RUN"])-1]
ind,=np.where((table["RUN"]==minrun))
mintime1= table['UNIXSTARTTIME'][ind]
mintime=mintime1.values
maxrun=table["RUN"][len(table["RUN"])-1]
now=time.time()
if (now-mintime/1000)<600:
minrun=table["RUN"][len(table["RUN"])-1]-1
print(now)
COL=[]
for i in DU:
DUx=[]
for j in AB:
for k in DOM:
try:
macaddress = db.doms.via_omkey((i,k), detid).dom_id
toas_all = sds.toashort(detid=detid, minrun=minrun, maxrun=maxrun, domid=macaddress, emitterid=j) # Prendere i dati basandosi sul tempo e non sul run????
QF_abdom=toas_all["QUALITYFACTOR"]
UTB_abdom=toas_all["UNIXTIMEBASE"]
TOAS_abdom=toas_all["TOA_S"]
UTB_abdom=UTB_abdom.values
up=np.where(UTB_abdom>(now-600))
down=np.where(UTB_abdom<(now))
intr=np.intersect1d(up,down)
UTB_abdom=UTB_abdom[intr]
QF_abdom=QF_abdom[intr]
QF_abdom=QF_abdom.values
QFlist=QF_abdom.tolist()
QFlist.sort(reverse=True)
QF_max=max(QF_abdom)
QF_max_index=np.where(QF_abdom==QF_max)
UTB_signal_min=UTB_abdom[QF_max_index]-80
UTB_signal_max=UTB_abdom[QF_max_index]+80
temp1=np.where(UTB_abdom>(UTB_signal_min[0]))
temp2=np.where(UTB_abdom<(UTB_signal_max[0]))
inter=np.intersect1d(temp1,temp2)
inter=inter.tolist()
signal_index=inter
QF_abdom_index=np.where(QF_abdom)
all_data_index=QF_abdom_index[0].tolist()
noise_index=diff(all_data_index,signal_index)
SIGNAL=QF_abdom[signal_index]
UTB_SIGNAL=UTB_abdom[signal_index]
NOISE=QF_abdom[noise_index]
NOISElist=NOISE.tolist()
NOISElist.sort(reverse=True)
noise_threshold=max(NOISE)
# First filter: 22 greatest
SIGNAL=SIGNAL.tolist()
SIGNAL_OLD=np.array(SIGNAL)
SIGNAL.sort(reverse=True)
QF_first=SIGNAL[0:22]
# Second filter: delete duplicates
QF_second=np.unique(QF_first)
QF_second=QF_second.tolist()
QF_second.sort(reverse=True)
# Third filter: If there are more than 11 elements I will eliminate the worst
if len(QF_second)>11:
QF_second=np.array(QF_second)
QF_third=[k for k in QF_second if (np.where(QF_second==k)[0][0]<11)]
else:
QF_third=QF_second
# Fourth filter: I remove the data if it is below the maximum noise
QF_fourth=[k for k in QF_third if k>(noise_threshold+(5*np.std(NOISE)))]
# Fifth filter: Check if the clicks are interspersed in the right way
# QF_fifth=[k for k in QF_fourth if (abs(k-max(QF_fourth))<abs(k-noise_threshold))]
Q=[]
for q in np.arange(len(QF_fourth)):
Q.append(np.where(SIGNAL_OLD==QF_fourth[q])[0][0])
UTB_fourth=np.array(UTB_SIGNAL.tolist())[Q]
UTB_fourth_l=UTB_fourth.tolist()
# UTB_fourth_l.sort()
D=[]
for g in np.arange(len(UTB_fourth_l)):
if ((np.mod((UTB_fourth_l[g]-UTB_fourth_l[0]),5)>0.5 and np.mod((UTB_fourth_l[g]-UTB_fourth_l[0]),5)<4.5) or (np.mod((UTB_fourth_l[g]-UTB_fourth_l[0]),5)>5)):
D.append(g)
for d in sorted(D, reverse=True):
del QF_fourth[d]
QF_fifth=QF_fourth
QF_OK=QF_fifth
# print(QF_OK)
NUM=len(QF_OK)
print(NUM)
if (NUM>7):
DUx.append(1.5)
elif (NUM<8 and NUM>3):
DUx.append(0.5)
elif (NUM<4 and NUM>0):
DUx.append(-0.5)
elif (NUM==0):
DUx.append(-1.5)
except:
stop=[]
DUx.append(-1.5)
COL.append(DUx)
fig = plt.figure()
ax = fig.add_subplot(111)
dom=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]
l=len(dom)
du1AB1=0.9*np.ones(l)
du1AB2=1*np.ones(l)
du1AB3=1.1*np.ones(l)
du2AB1=1.9*np.ones(l)
du2AB2=2*np.ones(l)
du2AB3=2.1*np.ones(l)
du3AB1=2.9*np.ones(l)
du3AB2=3*np.ones(l)
du3AB3=3.1*np.ones(l)
du4AB1=3.9*np.ones(l)
du4AB2=4*np.ones(l)
du4AB3=4.1*np.ones(l)
du5AB1=4.9*np.ones(l)
du5AB2=5*np.ones(l)
du5AB3=5.1*np.ones(l)
DU1=np.array(COL[0])
DU2=np.array(COL[1])
DU3=np.array(COL[2])
DU4=np.array(COL[3])
DU5=np.array(COL[4])
ind=np.where(DU2<1000)
iAB1=np.where(ind[0]<l)
iAB2_up=np.where(ind[0]>(l-1))
iAB2_down=np.where(ind[0]<2*l)
iAB2=np.intersect1d(iAB2_up,iAB2_down)
iAB3=np.where(ind[0]>(2*l-1))
colorsList = [(0, 0, 0),(1, 0.3, 0),(1, 1, 0),(0.2, 0.9, 0)]
CustomCmap = matplotlib.colors.ListedColormap(colorsList)
bounds=[-2,-1,0,1,2]
norma = colors.BoundaryNorm(bounds, CustomCmap.N)
color=ax.scatter(du1AB1,dom,s=20,c=DU1[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du1AB2,dom,s=20,c=DU1[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du1AB3,dom,s=20,c=DU1[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB1,dom,s=20,c=DU2[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB2,dom,s=20,c=DU2[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du2AB3,dom,s=20,c=DU2[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB1,dom,s=20,c=DU3[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB2,dom,s=20,c=DU3[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du3AB3,dom,s=20,c=DU3[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB1,dom,s=20,c=DU4[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB2,dom,s=20,c=DU4[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du4AB3,dom,s=20,c=DU4[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB1,dom,s=20,c=DU5[iAB1],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB2,dom,s=20,c=DU5[iAB2],norm=norma, marker = 's', cmap = CustomCmap);
color=ax.scatter(du5AB3,dom,s=20,c=DU5[iAB3],norm=norma, marker = 's', cmap = CustomCmap);
cbar = plt.colorbar(color)
cbar.ax.get_yaxis().set_ticks([])
for j, lab in enumerate(['$0. pings$','$1-3 pings$','$4-7 pings$','$>7. pings$']):
cbar.ax.text(3.5, (2 * j + 1) / 8.0, lab, ha='center', va='center')
cbar.ax.get_yaxis().labelpad = 18
matplotlib.pyplot.xticks(np.arange(1, 6, step=1))
matplotlib.pyplot.yticks(np.arange(0, 19, step=1))
matplotlib.pyplot.grid(color='k', linestyle='-', linewidth=0.2)
ax.set_xlabel('DUs', fontsize=18)
ax.set_ylabel('Floors', fontsize=18)
ts = now+3600
DATE=datetime.utcfromtimestamp(ts).strftime('%Y-%m-%d %H:%M:%S')
ax.set_title(r' %.16s Detection of the pings emitted by autonomous beacons' %DATE, fontsize=10)
# plt.show()
my_path = os.path.abspath(directory)
my_file = 'Online_Acoustic_Monitoring.png'
fig.savefig(os.path.join(my_path, my_file))
print(time.time())
check=False
check_time=time.time()-now
print(check_time)
time.sleep(abs(600-check_time))
check=True
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment