satnogs_waterfall_converter.py

satnogs_waterfall_converter.py

#!/usr/bin/env python3
import numpy as np
import h5py
import argparse
import os

if __name__ == "__main__":
    # Parse arguments
    parser = argparse.ArgumentParser(
        description="Convert a SatNOGS waterfall file to HDF5 format")
    parser.add_argument("-i", "--id", help="Observation ID", type=int)
    parser.add_argument("-t", "--starttime", help="Start time (YYYY-MM-DDTHH-MM-SS)")
    parser.add_argument("-f", "--freq", help="Center frequency (Hz)", type=float)
    parser.add_argument("-s", "--script_name", help="Script name")
    args = parser.parse_args()

    # Settings
    offset_in_stds = -2.0
    scale_in_stds = 8.0
    cfreq = args.freq
    obsid = args.id
    obs_timestamp = args.starttime
    script_name = args.script_name
    baud = 0
    tle = {"tle0": "", "tle1": "", "tle2": ""}
    fname = os.path.join(os.getenv("SATNOGS_OUTPUT_PATH"), "receiving_waterfall_%d_%s.dat" % (obsid, obs_timestamp))
    h5fname = os.path.join(os.getenv("SATNOGS_OUTPUT_PATH"), "observation_%d_%s.h5" % (obsid, obs_timestamp))

    # Read waterfall file
    fp = open(fname, "r")
    timestamp = np.fromfile(fp, dtype="|S32", count=1)[0]
    nchan = np.fromfile(fp, dtype='>i4', count=1)[0]
    samp_rate = np.fromfile(fp, dtype='>i4', count=1)[0]
    nfft_per_row = np.fromfile(fp, dtype='>i4', count=1)[0]
    center_freq = np.fromfile(fp, dtype='>f4', count=1)[0]
    endianness = np.fromfile(fp, dtype='>f4', count=1)[0] 


    data_dtypes = np.dtype([('tabs', 'int64'), ('spec', 'float32', (nchan, ))])
    data = np.fromfile(fp, dtype=data_dtypes)
    fp.close()

    nint = data.shape[0]

    tabs = data['tabs'] / 1000000.0
    trel = np.arange(nint) * nfft_per_row * nchan / float(samp_rate)

    waterfall = data['spec']
    
    freq = np.linspace(-0.5*samp_rate, 0.5*samp_rate, nchan, endpoint=False) / 1000.0

    # Compute time limits
    tmin, tmax = np.min(trel), np.max(trel)
    fmin, fmax = np.min(freq), np.max(freq)

    # Compute dynamic range limits
    c_idx = waterfall > -200.0
    if np.sum(c_idx) > 100:
        vmin = np.mean(waterfall[c_idx]) - 2.0 * np.std(waterfall[c_idx])
        vmax = np.mean(waterfall[c_idx]) + 6.0 * np.std(waterfall[c_idx])
    else:
        vmin = -100
        vmax = -50

    # Compute offset and scale
    waterfall_mean, waterfall_std = np.mean(waterfall, axis=0), np.std(waterfall, axis=0)
    waterfall_offset = waterfall_mean+offset_in_stds*waterfall_std
    waterfall_scale = scale_in_stds*waterfall_std/255.0

    # Convert waterfall to unsigned 8bit
    waterfall_8bit = np.clip((waterfall-waterfall_offset)/waterfall_scale, 0.0, 255.0).astype("uint8")

    # Create HDF5 file
    h5file = h5py.File(h5fname, "w")

    # Store observation attributes
    h5file.attrs["observation_id"] = obsid
    h5file.attrs["observation_timestamp"] = obs_timestamp
    h5file.attrs["center_frequency"] = cfreq
    h5file.attrs["center_frequency_unit"] = "Hz"
    h5file.attrs["script_name"] = script_name
    h5file.attrs["baud"] = baud
    h5file.attrs["tle_line_0"] = tle["tle0"]
    h5file.attrs["tle_line_1"] = tle["tle1"]
    h5file.attrs["tle_line_2"] = tle["tle2"]
    h5file.attrs["ground_station_id"] = int(os.getenv("SATNOGS_STATION_ID"))
    h5file.attrs["ground_station_lat"] = float(os.getenv("SATNOGS_STATION_LAT"))
    h5file.attrs["ground_station_lon"] = float(os.getenv("SATNOGS_STATION_LON"))
    h5file.attrs["ground_station_elev"] = float(os.getenv("SATNOGS_STATION_ELEV"))
    
    # Create waterfall group
    wf_group = h5file.create_group("waterfall")
    
    # Store waterfall attributes
    wf_group.attrs["start_time"] = timestamp
    wf_group.attrs["data_min"] = vmin
    wf_group.attrs["data_max"] = vmax
    wf_group.attrs["offset_in_stds"] = offset_in_stds
    wf_group.attrs["scale_in_stds"] = scale_in_stds
    
    # Store waterfall units
    wf_group.attrs["data_min_unit"] = "dB"
    wf_group.attrs["data_max_unit"] = "dB"
    wf_group.attrs["offset_unit"] = "dB"
    wf_group.attrs["scale_unit"] = "dB/div"
    wf_group.attrs["data_unit"] = "div"
    wf_group.attrs["relative_time_unit"] = "seconds"
    wf_group.attrs["absolute_time_unit"] = "seconds"
    wf_group.attrs["frequency_unit"] = "kHz"

    # Store waterfall datasets
    wf_group.create_dataset("offset", data=waterfall_offset, compression="gzip")
    wf_group.create_dataset("scale", data=waterfall_scale, compression="gzip")
    wf_group.create_dataset("data", data=waterfall_8bit, compression="gzip")
    wf_group.create_dataset("relative_time", data=trel, compression="gzip")
    wf_group.create_dataset("absolute_time", data=tabs, compression="gzip")
    wf_group.create_dataset("frequency", data=freq, compression="gzip")
    
    # Store waterfall labels
    wf_group["offset"].dims[0].label = "Time (seconds)"
    wf_group["scale"].dims[0].label = "Time (seconds)"
    wf_group["relative_time"].dims[0].label = "Time (seconds)"
    wf_group["absolute_time"].dims[0].label = "Time (seconds)"
    wf_group["frequency"].dims[0].label = "Frequency (kHz)"
    wf_group["data"].dims[0].label = "Frequency (kHz)"
    wf_group["data"].dims[1].label = "Time (seconds)"

    # Close file
    h5file.close()