2022-06-18 09:44:39 -06:00
#!/usr/bin/env python3
2022-06-18 10:15:35 -06:00
import sys
2022-06-21 02:49:28 -06:00
import argparse
2022-06-21 06:27:17 -06:00
import os
2022-06-26 10:40:11 -06:00
import re
2022-06-18 09:44:39 -06:00
import numpy as np
2022-06-21 10:03:38 -06:00
from strf . rfio import Spectrogram , get_site_info , get_frequency_info , get_satellite_info
2022-06-18 09:44:39 -06:00
import matplotlib . pyplot as plt
2022-06-18 10:01:38 -06:00
import matplotlib . dates as mdates
2022-06-18 10:24:45 -06:00
from matplotlib . backend_bases import MouseButton
2022-06-18 16:18:17 -06:00
from matplotlib . widgets import RectangleSelector
2022-06-18 16:37:37 -06:00
import matplotlib as mpl
2022-06-21 15:43:03 -06:00
from skyfield . api import EarthSatellite
from skyfield . api import load , wgs84 , utc
2022-06-22 16:33:05 -06:00
from datetime import datetime , timedelta
import h5py
import json
2022-06-26 10:40:11 -06:00
import requests
2022-07-31 20:18:00 -06:00
import shutil
2022-06-24 16:22:45 -06:00
from astropy . time import Time
2022-06-26 14:39:17 -06:00
from astropy . visualization import ZScaleInterval
2022-06-26 09:10:32 -06:00
C = 299792.458 # km/s
2022-06-20 12:01:21 -06:00
2022-06-20 11:36:03 -06:00
from modest import imshow
2022-06-19 09:25:11 -06:00
2022-06-22 06:46:39 -06:00
def main ( ) :
2022-06-25 02:08:23 -06:00
# Default settings
2022-06-25 09:27:18 -06:00
plt . style . use ( " dark_background " )
mpl . rcParams [ " keymap.save " ] . remove ( " s " )
mpl . rcParams [ " keymap.fullscreen " ] . remove ( " f " )
2022-06-22 16:33:05 -06:00
ts = load . timescale ( )
2022-06-21 02:49:28 -06:00
2022-06-25 02:08:23 -06:00
# Set defaults
if " ST_DATADIR " in os . environ :
site_fname = os . path . join ( os . environ [ " ST_DATADIR " ] , " data " , " sites.txt " )
freq_fname = os . path . join ( os . environ [ " ST_DATADIR " ] , " data " , " frequencies.txt " )
2022-06-26 10:31:45 -06:00
apikey_fname = os . path . join ( os . environ [ " ST_DATADIR " ] , " data " , " apikey.txt " )
2022-06-25 02:08:23 -06:00
else :
2022-06-26 10:42:12 -06:00
site_fname , freq_fname , apikey_fname = None , None , None
2022-06-25 02:08:23 -06:00
if " ST_TLEDIR " in os . environ :
tle_fname = os . path . join ( os . environ [ " ST_TLEDIR " ] , " bulk.tle " )
else :
tle_fname = None
# Parse input arguments
2022-06-25 01:21:29 -06:00
parser = argparse . ArgumentParser ( description = " rfplot: plot RF observations " , formatter_class = argparse . ArgumentDefaultsHelpFormatter )
2022-06-26 10:44:35 -06:00
parser . add_argument ( " -p " , " --path " , help = " File to read [STRF bin files as /path/to/file_??????.bin, /path/to/file_??????, /path/to/file, SatNOGS HDF5 artifact or number which indicates SatNOGs observation id] " , type = str , required = True )
2022-06-25 06:13:35 -06:00
parser . add_argument ( " -s " , " --start " , type = int , default = 0 , help = " Number of starting subintegration [STRF bin files] " )
parser . add_argument ( " -l " , " --length " , type = int , default = 3600 , help = " Number of subintegrations to plot [STRF bin files] " )
2022-07-31 19:26:14 -06:00
parser . add_argument ( " -C " , " --site " , type = int , help = " Site ID [STRF bin files] " , default = 9990 )
2022-06-25 02:08:23 -06:00
parser . add_argument ( " -F " , " --freqlist " , help = " List with frequencies " , default = freq_fname , type = str )
parser . add_argument ( " -d " , " --data " , help = " STRF dat path to load " , type = str )
parser . add_argument ( " -c " , " --catalog " , help = " TLE catalog to load " , default = tle_fname , type = str )
2022-06-21 02:49:28 -06:00
args = parser . parse_args ( )
2022-06-21 06:27:17 -06:00
2022-06-25 02:08:23 -06:00
# Input argument logic
if site_fname is None or not os . path . exists ( site_fname ) :
2022-06-21 06:49:10 -06:00
print ( f " Sites file not available under { site_fname } " )
2022-06-21 06:27:17 -06:00
sys . exit ( 1 )
2022-06-25 01:21:29 -06:00
site = get_site_info ( site_fname , args . site )
2022-06-21 06:49:10 -06:00
if site is None :
2022-06-25 01:21:29 -06:00
print ( f " Site with no: { args . site } does not exist " )
2022-06-21 06:27:17 -06:00
sys . exit ( 1 )
2022-06-25 02:08:23 -06:00
if not os . path . exists ( args . freqlist ) :
print ( f " Frequency list not available under { args . freqlist } " )
sys . exit ( 1 )
if not os . path . exists ( args . catalog ) :
print ( f " TLE catalog not available under { args . catalog } " )
2022-06-21 15:43:03 -06:00
2022-06-18 09:44:39 -06:00
# Read spectrogram
2022-06-26 10:31:45 -06:00
base_satellite = None
if re . match ( r " \ d+ " , args . path ) : # assume satnogs id
2022-06-26 10:42:12 -06:00
if apikey_fname is None :
print ( f " ST_DATADIR env variable not defined " )
sys . exit ( 1 )
2022-06-26 10:31:45 -06:00
if not os . path . exists ( apikey_fname ) :
print ( f " File containing API key not available under { apikey_fname } " )
sys . exit ( 1 )
with open ( apikey_fname , " r " ) as f :
apikey = f . read ( ) . strip ( )
config = requests . get ( f " https://db.satnogs.org/api/artifacts/?format=json&network_obs_id= { args . path } " , headers = { ' Authorization ' : f ' Token { apikey } ' } ) . json ( )
if len ( config ) == 0 :
print ( f " Observation with id: { args . path } has no artifact associated " )
sys . exit ( 1 )
filename = os . path . basename ( config [ 0 ] [ " artifact_file " ] )
args . path = filename
if not os . path . exists ( filename ) :
r = requests . get ( config [ 0 ] [ " artifact_file " ] , stream = True )
if r . status_code == 200 :
with open ( filename , ' wb ' ) as f :
r . raw . decode_content = True
shutil . copyfileobj ( r . raw , f )
elif r . status_code == 401 :
print ( f " Bad API key provided " )
sys . exit ( 1 )
2022-06-25 05:48:43 -06:00
fext = os . path . splitext ( args . path ) [ - 1 ]
if ( fext == " .h5 " ) or ( fext == " .hdf5 " ) :
s = Spectrogram . from_artifact ( args . path )
2022-06-22 17:24:22 -06:00
site = { " lat " : s . location [ " latitude " ] , " lon " : s . location [ " longitude " ] , " height " : s . location [ " altitude " ] }
site_location = wgs84 . latlon ( site [ " lat " ] , site [ " lon " ] , site [ " height " ] )
2022-06-26 06:05:33 -06:00
base_satellite = EarthSatellite ( s . tle [ - 2 ] , s . tle [ - 1 ] )
2022-06-22 17:24:22 -06:00
timestamps = [ x . replace ( tzinfo = utc ) for x in s . t ]
2022-06-26 06:05:33 -06:00
pos = ( base_satellite - site_location ) . at ( ts . utc ( timestamps ) )
2022-06-22 16:33:05 -06:00
_ , _ , _ , _ , _ , range_rate_base = pos . frame_latlon_and_rates ( site_location )
2022-06-22 17:24:22 -06:00
range_rate_base = range_rate_base . km_per_s
2022-06-25 05:48:43 -06:00
else :
s = Spectrogram . from_spectrogram ( args . path , args . start , args . length , args . site )
timestamps = [ x . replace ( tzinfo = utc ) for x in s . t ]
range_rate_base = [ 0 for x in timestamps ]
site_location = wgs84 . latlon ( site [ " lat " ] , site [ " lon " ] , site [ " height " ] )
2022-06-18 09:44:39 -06:00
# Create plot
2022-06-26 14:39:17 -06:00
vmin , vmax = ZScaleInterval ( ) . get_limits ( s . z )
2022-06-18 09:44:39 -06:00
2022-06-18 10:01:38 -06:00
# Time limits
tmin , tmax = mdates . date2num ( s . t [ 0 ] ) , mdates . date2num ( s . t [ - 1 ] )
# Frequency limits
2022-06-22 16:33:05 -06:00
fcen = s . fcen
2022-06-18 10:01:38 -06:00
fmin , fmax = ( s . freq [ 0 ] - fcen ) * 1e-6 , ( s . freq [ - 1 ] - fcen ) * 1e-6
2022-06-21 09:31:57 -06:00
frequencies = [ ]
2022-06-21 15:43:03 -06:00
satellite_info = [ ]
2022-06-25 05:48:43 -06:00
2022-06-25 02:08:23 -06:00
frequencies = get_frequency_info ( args . freqlist , fcen , s . freq [ 0 ] , s . freq [ - 1 ] )
2022-06-25 09:27:18 -06:00
names = ( " rise " , " culminate " , " set " )
2022-06-25 02:08:23 -06:00
t0 , t1 = ts . utc ( s . t [ 0 ] . replace ( tzinfo = utc ) ) , ts . utc ( s . t [ - 1 ] . replace ( tzinfo = utc ) )
satellite_info = get_satellite_info ( args . catalog , frequencies )
2022-06-21 10:03:38 -06:00
print ( f " Found { len ( frequencies ) } matching satellites " )
2022-06-21 09:31:57 -06:00
2022-06-18 16:18:17 -06:00
fig , ax = plt . subplots ( figsize = ( 10 , 6 ) )
2022-06-24 16:22:45 -06:00
2022-06-26 06:05:33 -06:00
def get_doppler_correction ( site , satellite , t , frequency ) : # frequency in Hz
_ , _ , _ , _ , _ , range_rate = ( satellite - site ) . at ( t ) . frame_latlon_and_rates ( site )
return range_rate . km_per_s / C * frequency
2022-06-24 16:22:45 -06:00
def plot_to_file ( array ) :
print ( array . shape )
ts1 = Time ( [ mdates . num2date ( x ) for x in array [ : , 0 ] ] ) . mjd
2022-06-26 10:38:35 -06:00
freqs = np . array ( [ x * 1e6 + fcen for x in array [ : , 1 ] ] )
2022-06-26 06:05:33 -06:00
if base_satellite is not None :
temp_t = ts . utc ( [ mdates . num2date ( x ) for x in array [ : , 0 ] ] )
freqs - = get_doppler_correction ( site_location , base_satellite , temp_t , fcen )
2022-06-24 16:22:45 -06:00
with open ( " mark.dat " , " w " ) as f :
for t , freq in zip ( ts1 , freqs ) :
2022-06-26 06:05:33 -06:00
f . write ( f " { t } { freq } 10 { args . site } \n " )
2022-06-24 16:22:45 -06:00
return ts1 , freqs
def file_to_plot ( ts1 , freqs ) :
2022-06-26 09:10:32 -06:00
if base_satellite is not None :
t = ts . utc ( [ x . replace ( tzinfo = utc ) for x in Time ( ts1 , format = ' mjd ' ) . datetime ] )
correction = get_doppler_correction ( site_location , base_satellite , t , fcen )
freqs + = correction
2022-06-24 16:22:45 -06:00
array = np . transpose ( np . array ( [
2022-06-25 09:27:18 -06:00
[ mdates . date2num ( x ) for x in Time ( ts1 , format = " mjd " ) . datetime ] ,
2022-06-24 16:22:45 -06:00
[ x - fcen * 1e-6 for x in freqs ]
] ) )
return array
2022-06-18 16:18:17 -06:00
mark = ax . scatter ( [ ] , [ ] , c = " white " , s = 5 )
2022-06-24 16:22:45 -06:00
2022-06-25 01:21:29 -06:00
if args . data is not None :
with open ( args . data , " r " ) as f :
2022-06-24 16:22:45 -06:00
lines = [ x . strip ( ) . split ( ) for x in f . readlines ( ) ]
mjds = [ float ( x [ 0 ] ) for x in lines ]
freqs = [ float ( x [ 1 ] ) for x in lines ]
array = file_to_plot ( mjds , freqs )
mark . set_offsets ( array )
2022-06-25 09:27:18 -06:00
line_fitting = ax . scatter ( [ ] , [ ] , edgecolors = " yellow " , s = 10 , facecolors = " none " )
2022-06-20 11:36:03 -06:00
# imshow(ax, s.z, vmin=vmin, vmax=vmax)
2022-06-21 15:43:03 -06:00
for sat_info in satellite_info :
satellite = EarthSatellite ( sat_info [ " tle " ] [ - 2 ] , sat_info [ " tle " ] [ - 1 ] )
2022-06-21 16:55:52 -06:00
t , events = satellite . find_events ( site_location , t0 , t1 , altitude_degrees = 0.0 )
2022-06-25 09:02:03 -06:00
alt , _ , _ = ( satellite - site_location ) . at ( t0 ) . altaz ( )
if len ( t ) == 0 and alt . degrees > 0 :
t , events = [ t0 , t1 ] , [ 0 , 2 ]
2022-06-21 15:43:03 -06:00
if len ( t ) > 0 :
pairs = [ ( ti , event ) for ti , event in zip ( t , events ) ]
if pairs [ 0 ] [ 1 ] in [ 1 , 2 ] :
pairs = [ ( t0 , 0 ) ] + pairs # pad with rise
if pairs [ - 1 ] [ 1 ] in [ 0 , 1 ] :
pairs = pairs + [ ( t1 , 2 ) ] # pad with set
pairs = [ ( ti , event ) for ti , event in pairs if event != 1 ] # remove culminations
2022-06-21 16:46:22 -06:00
sat_info [ " timeslot " ] = [ ( pairs [ i ] [ 0 ] . utc_datetime ( ) , pairs [ i + 1 ] [ 0 ] . utc_datetime ( ) ) for i in range ( 0 , len ( pairs ) , 2 ) ]
2022-06-21 15:43:03 -06:00
for timeslot in sat_info [ " timeslot " ] :
2022-06-22 17:24:22 -06:00
selected_pairs = [ x for x in zip ( timestamps , range_rate_base ) if x [ 0 ] > = timeslot [ 0 ] and x [ 0 ] < = timeslot [ 1 ] ]
2022-06-22 16:33:05 -06:00
selected_timestamps = [ x [ 0 ] for x in selected_pairs ]
selected_range_rate_base = np . array ( [ x [ 1 ] for x in selected_pairs ] )
2022-06-21 16:46:22 -06:00
pos = ( satellite - site_location ) . at ( ts . utc ( selected_timestamps ) )
_ , _ , _ , _ , _ , range_rate = pos . frame_latlon_and_rates ( site_location )
2022-06-22 16:33:05 -06:00
range_rate_sat = range_rate . km_per_s
2022-06-21 16:46:22 -06:00
2022-06-22 16:33:05 -06:00
for fsat in sat_info [ " frequencies " ] :
fbase = fcen * 1e-6
dfreq = ( 1 - range_rate_sat / C ) * fsat - ( 1 - selected_range_rate_base / C ) * fbase
2022-06-22 03:22:12 -06:00
tt = [ mdates . date2num ( x ) for x in selected_timestamps ]
2022-06-22 16:33:05 -06:00
ax . plot ( tt , dfreq , c = " orange " )
ax . text ( tt [ 0 ] , dfreq [ 0 ] , sat_info [ " noradid " ] , c = " orange " )
2022-06-21 15:43:03 -06:00
2022-06-21 04:07:06 -06:00
image = imshow ( ax , s . z , origin = " lower " , aspect = " auto " , interpolation = " None " ,
2022-06-18 10:01:38 -06:00
vmin = vmin , vmax = vmax ,
extent = [ tmin , tmax , fmin , fmax ] )
2022-06-22 03:22:12 -06:00
2022-06-20 11:36:03 -06:00
mode = {
2022-06-21 04:07:06 -06:00
" current_mode " : None ,
" vmin " : vmin ,
" vmax " : vmax
2022-06-20 11:36:03 -06:00
}
2022-06-21 04:07:06 -06:00
2022-06-18 16:18:17 -06:00
def line_select_callback ( eclick , erelease ) :
x1 , y1 = eclick . xdata , eclick . ydata
x2 , y2 = erelease . xdata , erelease . ydata
2022-06-20 11:36:03 -06:00
if mode [ " current_mode " ] == " fit " :
t1_ind = round ( len ( s . t ) * ( x1 - tmin ) / ( tmax - tmin ) )
t2_ind = round ( len ( s . t ) * ( x2 - tmin ) / ( tmax - tmin ) )
f1_ind = round ( len ( s . freq ) * ( y1 - fmin ) / ( fmax - fmin ) )
f2_ind = round ( len ( s . freq ) * ( y2 - fmin ) / ( fmax - fmin ) )
2022-06-21 02:49:28 -06:00
submat = s . z [ f1_ind : f2_ind , t1_ind : t2_ind ]
# TODO perform some action on submat
2022-06-20 11:36:03 -06:00
elif mode [ " current_mode " ] == " delete " :
array = mark . get_offsets ( )
maskx = np . logical_and ( array [ : , 0 ] > = min ( x1 , x2 ) , array [ : , 0 ] < = max ( x1 , x2 ) )
masky = np . logical_and ( array [ : , 1 ] > = min ( y1 , y2 ) , array [ : , 1 ] < = max ( y1 , y2 ) )
mask = np . logical_and ( maskx , masky )
mark . set_offsets ( array [ np . logical_not ( mask ) , : ] )
2022-06-18 16:18:17 -06:00
fig . canvas . draw ( )
2022-06-21 02:53:05 -06:00
current_mode = mode [ " current_mode " ]
print ( f " select over { x1 } , { y1 } , { x2 } , { y2 } in { current_mode } mode " )
2022-06-18 16:18:17 -06:00
2022-06-25 09:27:18 -06:00
selector = RectangleSelector ( ax , line_select_callback , useblit = True , button = [ 1 ] , minspanx = 5 , minspany = 5 , spancoords = " pixels " , props = { " edgecolor " : " white " , " fill " : False } )
2022-06-18 16:18:17 -06:00
selector . active = False
2022-06-20 11:36:03 -06:00
2022-06-18 10:01:38 -06:00
ax . xaxis_date ( )
date_format = mdates . DateFormatter ( " %F \n % H: % M: % S " )
ax . xaxis . set_major_formatter ( date_format )
fig . autofmt_xdate ( rotation = 0 , ha = " center " )
ax . set_xlabel ( " Time (UTC) " )
ax . set_ylabel ( f " Frequency (MHz) - { fcen * 1e-6 : g } MHz " )
2022-06-18 16:18:17 -06:00
def add_point ( scatter , point ) :
array = scatter . get_offsets ( )
array = np . vstack ( [ array , point ] )
2022-06-24 16:22:45 -06:00
plot_to_file ( array )
2022-06-18 16:18:17 -06:00
scatter . set_offsets ( array )
fig . canvas . draw ( )
def handle ( key , x , y ) :
print ( f " pressed { key } over x= { x } y= { y } " )
if key == " d " :
selector . active = True
2022-06-20 11:36:03 -06:00
mode [ " current_mode " ] = " delete "
2022-06-18 16:37:37 -06:00
elif key == " s " :
point = ( x , y )
add_point ( line_fitting , point )
elif key == " f " :
print ( " performing fitting on " )
2022-06-20 11:36:03 -06:00
mode [ " current_mode " ] = " fit "
selector . active = True
2022-06-18 16:37:37 -06:00
print ( line_fitting . get_offsets ( ) )
2022-06-18 16:56:54 -06:00
elif key == " r " :
print ( " performing reset " )
2022-06-20 11:36:03 -06:00
mode [ " current_mode " ] = None
selector . active = False
2022-06-18 16:56:54 -06:00
mark . set_offsets ( np . empty ( ( 0 , 2 ) , float ) )
line_fitting . set_offsets ( np . empty ( ( 0 , 2 ) , float ) )
fig . canvas . draw ( )
2022-06-21 04:07:06 -06:00
elif key == " b " :
mode [ " vmax " ] * = 0.95
image . set_clim ( mode [ " vmin " ] , mode [ " vmax " ] )
fig . canvas . draw ( )
elif key == " v " :
mode [ " vmax " ] * = 1.05
image . set_clim ( mode [ " vmin " ] , mode [ " vmax " ] )
fig . canvas . draw ( )
2022-06-18 16:37:37 -06:00
sys . stdout . flush ( )
2022-06-18 16:18:17 -06:00
def on_press ( event ) :
handle ( event . key , event . xdata , event . ydata )
sys . stdout . flush ( )
def on_click ( event ) :
if event . button is MouseButton . MIDDLE :
point = ( event . xdata , event . ydata )
add_point ( mark , point )
print ( f " { event . xdata } { fcen + event . ydata } " )
sys . stdout . flush ( )
2022-06-25 09:27:18 -06:00
fig . canvas . mpl_connect ( " key_press_event " , on_press )
fig . canvas . mpl_connect ( " button_press_event " , on_click )
2022-06-18 09:44:39 -06:00
plt . show ( )
2022-06-22 06:46:39 -06:00
if __name__ == " __main__ " :
main ( )