From 78b99dda5a3be8fb61d938c478a79bd56ea90440 Mon Sep 17 00:00:00 2001 From: Cees Bassa Date: Sat, 25 Jun 2022 13:48:43 +0200 Subject: [PATCH] Move artifact definition to Spectrogram class --- rfplot.py | 35 +++++++------------------- strf/rfio.py | 71 ++++++++++++++++++++++++++++++++++++++-------------- 2 files changed, 61 insertions(+), 45 deletions(-) diff --git a/rfplot.py b/rfplot.py index fc80605..8fce783 100755 --- a/rfplot.py +++ b/rfplot.py @@ -20,25 +20,6 @@ from astropy.time import Time from modest import imshow -class Artifact: - def __init__(self, filename): - hdf5_file = h5py.File(filename, 'r') - - if hdf5_file.attrs['artifact_version'] != 2: - raise Exception(hdf5_file.attrs['artifact_version']) - - wf = hdf5_file.get('waterfall') - metadata = json.loads(hdf5_file.attrs['metadata']) - start_time = datetime.strptime(wf.attrs['start_time'].decode('ascii'), '%Y-%m-%dT%H:%M:%S.%fZ') - print(metadata) - self.z = np.transpose((np.array(wf['data']) * np.array(wf['scale']) + np.array(wf['offset']))) - self.fcen = float(metadata['frequency']) - self.freq = np.array(hdf5_file.get('waterfall').get("frequency")) + self.fcen - self.t = [ start_time + timedelta(seconds=x) for x in np.array(hdf5_file.get('waterfall').get("relative_time"))] - self.location = metadata["location"] - self.tle = [ x for x in metadata["tle"].split("\n") if x.strip() != "" ] - - def main(): # Default settings plt.style.use('dark_background') @@ -84,13 +65,9 @@ def main(): print(f"TLE catalog not available under {args.catalog}") # Read spectrogram - if args.artifact is None: - s = 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"]) - else: - s = Artifact(args.artifact) + fext = os.path.splitext(args.path)[-1] + if (fext == ".h5") or (fext == ".hdf5"): + s = Spectrogram.from_artifact(args.path) site = {"lat" : s.location["latitude"],"lon" : s.location["longitude"],"height" : s.location["altitude"]} site_location = wgs84.latlon(site["lat"], site["lon"], site["height"]) satellite = EarthSatellite(s.tle[-2], s.tle[-1]) @@ -98,6 +75,11 @@ def main(): pos = (satellite - site_location).at(ts.utc(timestamps)) _, _, _, _, _, range_rate_base = pos.frame_latlon_and_rates(site_location) range_rate_base = range_rate_base.km_per_s + 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"]) # Create plot vmin, vmax = np.percentile(s.z, (5, 99.95)) @@ -111,6 +93,7 @@ def main(): frequencies = [] satellite_info = [] + frequencies = get_frequency_info(args.freqlist, fcen, s.freq[0], s.freq[-1]) names = ('rise', 'culminate', 'set') t0,t1 = ts.utc(s.t[0].replace(tzinfo=utc)), ts.utc(s.t[-1].replace(tzinfo=utc)) diff --git a/strf/rfio.py b/strf/rfio.py index b64943b..9fc3e55 100644 --- a/strf/rfio.py +++ b/strf/rfio.py @@ -4,6 +4,9 @@ import re import sys import glob +import h5py +import json + import numpy as np from datetime import datetime, timedelta @@ -11,7 +14,20 @@ from datetime import datetime, timedelta class Spectrogram: """Spectrogram class""" - def __init__(self, froot, ifile, nsub, siteid): + def __init__(self, t, fcen, bw, freq, z, siteid, location, tle): + self.z = z + self.nchan = self.z.shape[0] + self.nsub = self.z.shape[1] + self.t = t + self.fcen = fcen + self.bw = bw + self.freq = freq + self.siteid = siteid + self.location = location + self.tle = tle + + @classmethod + def from_spectrogram(cls, froot, ifile, nsub, siteid): """Define a spectrogram""" path, prefix = extract_path_and_prefix(froot) @@ -67,17 +83,35 @@ class Spectrogram: isub += 1 ifile += 1 - self.z = np.transpose(np.vstack(zs)) - self.t = t - self.freq = freq - self.fcen = freq_cen - self.bw = bw - self.siteid = siteid - self.nchan = self.z.shape[0] - self.nsub = self.z.shape[1] + z = np.transpose(np.vstack(zs)) + + print(f"Read spectrogram\n{nchan} channels, {nsub} subints\nFrequency: {freq_cen * 1e-6:g} MHz\nBandwidth: {bw * 1e-6:g} MHz") - print(f"Read spectrogram\n{self.nchan} channels, {self.nsub} subints\nFrequency: {self.fcen * 1e-6:g} MHz\nBandwidth: {self.bw * 1e-6:g} MHz") + return cls(t, freq_cen, bw, freq, z, siteid, None, None) + + @classmethod + def from_artifact(cls, filename): + hdf5_file = h5py.File(filename, "r") + + if hdf5_file.attrs["artifact_version"] != 2: + raise Exception(hdf5_file.attrs["artifact_version"]) + + wf = hdf5_file.get("waterfall") + metadata = json.loads(hdf5_file.attrs["metadata"]) + start_time = datetime.strptime(wf.attrs["start_time"].decode("ascii"), "%Y-%m-%dT%H:%M:%S.%fZ") + z = np.transpose((np.array(wf["data"]) * np.array(wf["scale"]) + np.array(wf["offset"]))) + fcen = float(metadata["frequency"]) + freq = np.array(hdf5_file.get("waterfall").get("frequency")) + fcen + bw = (freq[1] - freq[0]) * z.shape[0] + t = [ start_time + timedelta(seconds=x) for x in np.array(hdf5_file.get("waterfall").get("relative_time"))] + location = metadata["location"] + tle = [ x for x in metadata["tle"].split("\n") if x.strip() != "" ] + nsub, nchan = z.shape + + print(f"Read artifact\n{nchan} channels, {nsub} subints\nFrequency: {fcen * 1e-6:g} MHz\nBandwidth: {bw * 1e-6:g} MHz") + + return cls(t, fcen, bw, freq, z, None, location, tle) def extract_path_and_prefix(fname): basename, dirname = os.path.basename(fname), os.path.dirname(fname) @@ -94,21 +128,21 @@ def extract_path_and_prefix(fname): return dirname, prefix def parse_header(header_b): - header_s = header_b.decode('ASCII').strip('\x00') + header_s = header_b.decode("ASCII").strip("\x00") regex = r"^HEADER\nUTC_START (.*)\nFREQ (.*) Hz\nBW (.*) Hz\nLENGTH (.*) s\nNCHAN (.*)\nNSUB (.*)\nEND\n$" try: match = re.fullmatch(regex, header_s, re.MULTILINE) except AttributeError: match = re.match(regex, header_s, flags=re.MULTILINE) - utc_start = datetime.strptime(match.group(1), '%Y-%m-%dT%H:%M:%S.%f') + utc_start = datetime.strptime(match.group(1), "%Y-%m-%dT%H:%M:%S.%f") - return {'utc_start': utc_start, - 'freq': float(match.group(2)), - 'bw': float(match.group(3)), - 'length': float(match.group(4)), - 'nchan': int(match.group(5)), - 'nsub': int(match.group(6))} + return {"utc_start": utc_start, + "freq": float(match.group(2)), + "bw": float(match.group(3)), + "length": float(match.group(4)), + "nchan": int(match.group(5)), + "nsub": int(match.group(6))} def get_site_info(fname, site_id): @@ -143,7 +177,6 @@ def get_frequency_info(fname, fcen, fmin, fmax): C = 299792.458 # km/s padding = fcen * 20 / C # padding in MHz fmin_padding,fmax_padding = (fmin - padding) * 1e-6, (fmax + padding) * 1e-6 - frequencies = {} for line in lines: if "#" in line: