From 098406bdd60e52d6a69446a9a802e98e9439f928 Mon Sep 17 00:00:00 2001 From: "Fabian P. Schmidt" Date: Sat, 20 Nov 2021 06:38:23 +0100 Subject: [PATCH] Add SatNOGS Artifact Analysis (WIP) Signed-off-by: Fabian P. Schmidt --- contrib/README.md | 41 ++++++- contrib/analyze_artifact.py | 158 +++++++++++++++++++++++++++ contrib/download_satnogs_artifact.py | 18 ++- contrib/download_satnogs_tle.py | 26 +++++ contrib/requirements.txt | 4 +- contrib/test_artifact_download.sh | 7 ++ 6 files changed, 242 insertions(+), 12 deletions(-) create mode 100755 contrib/analyze_artifact.py create mode 100755 contrib/download_satnogs_tle.py create mode 100755 contrib/test_artifact_download.sh diff --git a/contrib/README.md b/contrib/README.md index c502e0e..f29a5cb 100644 --- a/contrib/README.md +++ b/contrib/README.md @@ -20,14 +20,43 @@ $ ./contrib/find_good_satnogs_artifacts.py ``` ## Download SatNOGS Artifacts + ``` -$ ./contrib/download_satnogs_artifact.py 2786575 -Artifact Metadata for Observation #2786575 found. -Download failed for https://db-satnogs.freetls.fastly.net/media/artifacts/b4975058-04eb-4ab7-9c40-9bcce76d94db.h5 +$ ./contrib/download_satnogs_artifact.py 4950356 +Artifact Metadata for Observation #4950356 found. +Artifact saved in /home/pi/data/artifacts/4950356.h5 +``` + +## Download SatNOGS Observation TLEs + +To add the TLE used in a specific SatNOGS Observation to your catalog, the following command can be used: +``` +$ ./contrib/download_satnogs_tle.py 4950356 +TLE saved in /home/pi/data/tles/satnogs/4950356.txt +``` + +## Plot SatNOGS Artifacts + +This can be done using [spectranalysis](https://github.com/kerel-fs/spectranalysis). + +## Ultra-Short Analysis Guide +``` +OBSERVATION_ID=4950356 +./contrib/download_satnogs_tle.py $OBSERVATION_ID +./contrib/download_satnogs_artifact.py $OBSERVATION_ID +./contrib/analyze_artifact.py --observation_id $OBSERVATION_ID --site_id 977 +rffit -d "$SATNOGS_DOPPLER_OBS_DIR/$OBSERVATION_ID.dat" -c "$SATNOGS_TLE_DIR/$OBSERVATION_ID.txt" -i 27844 -s 977 ``` ``` -$ ./contrib/download_satnogs_artifact.py 4443137 -Artifact Metadata for Observation #4443137 found. -Artifact for Observation #4443137 saved in '/tmp/tmp1rdpnz_k' +$ ./contrib/download_satnogs_tle.py 4950356 +TLE saved in /mnt/old_home/kerel/c4/satnogs/data/tles/satnogs/4950356.txt +$ ./contrib/download_satnogs_artifact.py 4950356 +Artifact Metadata for Observation #4950356 found. +Artifact saved in /mnt/old_home/kerel/c4/satnogs/data/artifacts/4950356.h5 +$ ./contrib/analyze_artifact.py --observation_id 4950356 --site_id 977 +Load /mnt/old_home/kerel/c4/satnogs/data/artifacts/4950356.h5 +Extract measurements... +Data written in /mnt/old_home/kerel/c4/satnogs/data/doppler_obs/4950356.dat +$ rffit -d ${SATNOGS_DOPPLER_OBS_DIR}/4950356.dat -c ${SATNOGS_TLE_DIR}/4950356.txt -i 27844 -s 977 ``` diff --git a/contrib/analyze_artifact.py b/contrib/analyze_artifact.py new file mode 100755 index 0000000..cce8b20 --- /dev/null +++ b/contrib/analyze_artifact.py @@ -0,0 +1,158 @@ +#!/usr/bin/env python3 + +from astropy.time import Time + +from datetime import datetime, timedelta + +import argparse +import ephem +import h5py +import json +import os + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + + +def load_artifact(filename): + hdf5_file = h5py.File(filename, 'r') + if hdf5_file.attrs['artifact_version'] != 2: + print("unsupported artifact version {}".format(hdf5_file.attrs['artifact_version'])) + # return + wf = hdf5_file.get('waterfall') + data = (np.array(wf['data']) * np.array(wf['scale']) + np.array(wf['offset'])) + metadata = json.loads(hdf5_file.attrs['metadata']) + + return hdf5_file, wf, data, metadata + +def extract_peaks(data): + """ + Returns + ------- + points: (time_index, freq_index) + measurements: measurement tuples of relative time in seconds and relative frequency in hertz + """ + + snr = (np.max(data, axis=1) - np.mean(data, axis=1)) / np.std(data, axis=1) + + # Integrate each channel along the whole observation, + # This filter helps under the assumption of minimal Doppler deviaion + snr_integrated = np.zeros(data[0].shape) + for row in data: + snr_integrated += (row - np.mean(row)) / np.std(row) + snr_integrated = pd.Series(snr_integrated/len(snr_integrated)) + + SNR_CUTOFF = 4 + CHANNEL_WINDOW_SIZE = 4 + channel_mask = (snr_integrated.diff().abs() / snr_integrated.diff().std() > SNR_CUTOFF).rolling(CHANNEL_WINDOW_SIZE, min_periods=1).max() + rows=[] + for i, row in enumerate(channel_mask): + if not row: + continue + rows.append(i) + + # Select only maximum values in masked channels + points = [] + for i,x in enumerate(np.argmax(data, axis=1)): + if not x in rows: + continue + points.append([i, x]) + points = np.array(points) + measurements = np.vstack((wf['relative_time'][points[:,0]], + [wf['frequency'][x] for x in points[:,1]])) + return snr, measurements, points + +def plot_measurements(wf, data, measurements, snr): + plt.plot(wf['relative_time'][:], + [wf['frequency'][x] for x in np.argmax(data[:], axis=1)], + '.', + label='all') + plt.plot(measurements[:,0], + measurements[:,1], + '.', + label='automatic channel mask') + plt.plot(wf['relative_time'][2.4 < snr], + [wf['frequency'][x] for x in np.argmax(data[2.4 < snr], axis=1)], + '.', + label="SNR > 2.4") + plt.legend() + plt.title("Observation #4991792 - Maximum Values") + plt.grid() + plt.xlabel('Elapsed Time / s') + plt.ylabel('rel. Frequency / kHz') + plt.show() + + + +def dedoppler(measurements, metadata): + tle = metadata['tle'].split('\n') + start_time = datetime.strptime(wf.attrs['start_time'].decode('ascii'), + '%Y-%m-%dT%H:%M:%S.%fZ') + f_center = float(metadata['frequency']) + + # Initialize SGP4 propagator / pyephem + satellite = ephem.readtle('sat', tle[1], tle[2]) + observer = ephem.Observer() + observer.lat = str(metadata['location']['latitude']) + observer.lon = str(metadata['location']['longitude']) + observer.elevation = metadata['location']['altitude'] + + def remove_doppler_correction(t, freq): + """ + Arguments + --------- + t - float: Time in seconds + freq - float: Relative Frequency in herz + """ + observer.date = t + satellite.compute(observer) + v = satellite.range_velocity + df = f_center * v / ephem.c + return f_center + freq - df*2 + + output = [] + for dt,df in measurements.T: + t = start_time + timedelta(seconds=dt) + f = f_center + df + freq_recv = remove_doppler_correction(t, f) + output.append((t, freq_recv)) + return output + +def save_rffit_data(filename, measurements, site_id): + with open(filename, 'w') as file_out: + for time, freq in measurements: + line = '{:.6f}\t{:.2f}\t1.0\t{}\n'.format(Time(time).mjd, freq, site_id) + file_out.write(line) + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Analyze SatNOGS Artifact.') + parser.add_argument('--observation_id', type=str, + help='SatNOGS Observation ID') + parser.add_argument('--filename', type=str, + help='SatNOGS Artifacts File') + parser.add_argument('--output_file', type=str, + help='STRF-compatible output file') + parser.add_argument('--site_id', type=int, required=True, + help='STRF Site ID') + args = parser.parse_args() + + if args.observation_id: + # Use canonic file paths, ignoring `--filename` and `--output_path` + filename = '{}/{}.h5'.format(os.getenv('SATNOGS_ARTIFACTS_DIR'), args.observation_id) + output_file = '{}/{}.dat'.format(os.getenv('SATNOGS_DOPPLER_OBS_DIR'), args.observation_id) + else: + if not any([args.filename, args.output_file]): + print('ERROR: Missing arguments') + filename = args.filename + output_file = args.output_file + + print('Load {}'.format(filename)) + hdf5_file, wf, data, metadata = load_artifact(filename) + + print('Extract measurements...') + snr, measurements, points = extract_peaks(data) + # plot_measurements(wf, data, measurements, snr) + m2 = dedoppler(measurements, metadata) + save_rffit_data(output_file, m2, site_id=args.site_id) + print('Data written in {}'.format(output_file)) diff --git a/contrib/download_satnogs_artifact.py b/contrib/download_satnogs_artifact.py index 50d4447..81c9866 100755 --- a/contrib/download_satnogs_artifact.py +++ b/contrib/download_satnogs_artifact.py @@ -1,11 +1,13 @@ #!/usr/bin/env python3 import argparse -import tempfile -import sys import logging import requests import settings +import shutil +import sys +import tempfile +import os from urllib.parse import urljoin from pprint import pprint @@ -41,7 +43,7 @@ def fetch_artifact(url, artifact_filename): fname.write(chunk) -def download_artifact(observation_id): +def download_artifact(observation_id, filename): try: artifact_metadata = fetch_artifact_metadata(network_obs_id=observation_id) except requests.HTTPError: @@ -57,8 +59,13 @@ def download_artifact(observation_id): try: artifact_file_url = artifact_metadata[0]['artifact_file'] artifact_file = tempfile.NamedTemporaryFile(delete=False) + fetch_artifact(artifact_file_url, artifact_file.name) - print("Artifact for Observation #{} saved in '{}'".format(observation_id, artifact_file.name)) + + artifact_file.close() + shutil.copy(artifact_file.name, filename) + os.remove(artifact_file.name) + print("Artifact saved in {}".format(filename)) except requests.HTTPError: print('Download failed for {}'.format(artifact_file_url)) return @@ -74,4 +81,5 @@ if __name__ == "__main__": logging.basicConfig(level=logging.INFO) for observation_id in args.observation_ids: - download_artifact(observation_id) + download_artifact(observation_id, + '{}/{}.h5'.format(os.getenv('SATNOGS_ARTIFACTS_DIR'), observation_id)) diff --git a/contrib/download_satnogs_tle.py b/contrib/download_satnogs_tle.py new file mode 100755 index 0000000..8cc460a --- /dev/null +++ b/contrib/download_satnogs_tle.py @@ -0,0 +1,26 @@ +#!/usr/bin/env python3 + +import argparse +import os + +from satnogs_api_client import fetch_observation_data + + +if __name__ == '__main__': + parser = argparse.ArgumentParser(description='Download TLE from a specific Observation in SatNOGS Network.') + parser.add_argument('observation_id', type=int, + help='SatNOGS Observation ID') + args = parser.parse_args() + + obs = fetch_observation_data([args.observation_id])[0] + + filename = '{}/{}.txt'.format(os.getenv('SATNOGS_TLE_DIR'), args.observation_id) + + with open(filename, 'w') as f: + f.write(obs['tle0']) + f.write('\n') + f.write(obs['tle1']) + f.write('\n') + f.write(obs['tle2']) + f.write('\n') + print("TLE saved in {}".format(filename)) diff --git a/contrib/requirements.txt b/contrib/requirements.txt index 358f2c7..eda36ad 100644 --- a/contrib/requirements.txt +++ b/contrib/requirements.txt @@ -1,8 +1,10 @@ astropy -ephem matplotlib numpy Pillow python-decouple requests git+https://gitlab.com/librespacefoundation/satnogs/python-satnogs-api.git@e20a7d3c +h5py~=3.6.0 +pandas~=1.3.4 +ephem~=4.3.1 diff --git a/contrib/test_artifact_download.sh b/contrib/test_artifact_download.sh new file mode 100755 index 0000000..395253f --- /dev/null +++ b/contrib/test_artifact_download.sh @@ -0,0 +1,7 @@ +#!/usr/bin/bash + +DB_API_TOKEN="4f20a493a3f5fd85074d61db944365bf94987613" +URL="https://db-satnogs.freetls.fastly.net/media/artifacts/b4975058-04eb-4ab7-9c40-9bcce76d94db.h5" +echo "Authorization: Token $DB_API_TOKEN" +curl -H "Authorization: Token $DB_API_TOKEN" "$URL" +