Unverified Commit 6bc759b4 authored by Tammo Jan Dijkema's avatar Tammo Jan Dijkema
Browse files

Add preliminary map view

parent e3ebfaf1
#!/usr/bin/env python3
"""
Receive numbers from gnuradio (over UDP) at about 2 per second, plot them
on a map
Tammo Jan Dijkema, July 2021
"""
import time
import matplotlib.pyplot as plt
import pandas as pd
import zmq
import re
import struct
import numpy as np
from telescope import Telescope
from astropy.coordinates import get_sun, SkyCoord, AltAz, EarthLocation, Galactic
from astropy.time import Time
import astropy.units as u
from matplotlib import cm
from threading import Thread
plt.ion()
import logging
logger = logging.getLogger()
logger.setLevel(level=logging.INFO)
stoprequest = False
def set_stoprequest(event):
global stoprequest
stoprequest = True
signal = 0.
prev_signal = 0.
datafile = open(f"data{int(Time.now().unix)}.txt", "w")
print("Time (UTC),separation (deg),azimuth from north (deg), elevation (deg), signal, source", file=datafile)
def signal_thread():
global signal
context = zmq.Context()
receiver = context.socket(zmq.PULL)
receiver.connect("tcp://127.0.0.1:5557")
while True:
buff = receiver.recv()
signal = np.frombuffer(buff, dtype="float32")[-1]
s = Thread(target=signal_thread)
s.daemon = True
s.start()
dt = Telescope(consoleHost='localhost')
data = np.zeros((180, 360), dtype=float)
npoints = np.zeros(data.shape, dtype='i8')
# create plot
plt.ion()
fig, ax = plt.subplots()
fig.canvas.mpl_connect('close_event', set_stoprequest)
plot = ax.imshow(data, origin='lower', extent=(-180, 180, -90, 90))
ax.set_xlabel("Galactische longitude (graden)")
ax.set_ylabel("Galactische latitude (graden)")
plt.show()
plt.pause(0.001)
prev_points = None
maxsignal = 0.
minsignal = 1.0e9
while not stoprequest:
plt.pause(0.2)
if signal == prev_signal:
continue
prev_signal = signal
pos = dt.radec
if pos is None:
print("Telescope not ready yet")
continue
pos = pos.transform_to(Galactic)
l = round(pos.l.deg)
b = round(pos.b.deg)
if prev_points is not None:
try:
prev_points.remove()
except:
pass
#print(Time.now(), separation, az.to(u.deg).value, alt.to(u.deg).value, signal, clean_name(closest_source), file=datafile, flush=True, sep=',')
if signal > maxsignal:
maxsignal = signal
if signal < minsignal:
minsignal = signal
n = npoints[b + 90, l + 180]
data[b + 90, l + 180] = (signal + n * data[b + 90, l + 180]) / (n + 1)
npoints[b + 90, l + 180] += 1
plot.set_data(data)
plot.set_clim(vmin=minsignal, vmax=maxsignal)
fig.canvas.flush_events()
xmin, xmax = ax.get_xlim()
ymin, ymax = ax.get_ylim()
prev_points = ax.scatter([l+0.5], [b+0.5], s=100, facecolors='none', edgecolors='r')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
# show the plot
plt.show()
plt.pause(0.0001)
Markdown is supported
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