Skip to content
Snippets Groups Projects

Add acoustics monitoring

Merged Add acoustics monitoring
cguidi/km3mon:undefined into master
Merged Carlo Guidi requested to merge cguidi/km3mon:undefined into master
Compare
and
1 file
+ 273
− 0
Compare changes
  • Side-by-side
  • Inline
scripts/acoustics.py 0 → 100644
+ 273
− 0
#!/usr/bin/env python
# coding=utf-8
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