read Madrigal GNSS LOS file

read_gnss_los.py

"""
read Madrigal GNSS LOS file
based on http://cedar.openmadrigal.org/static/siteSpecific/programming_tips.pdf
"""
import time
import h5py
import numpy
import argparse


def getTimeData(madrigalFile: str, timeArr, unixTime):
    """getTimeData returns a numpy recarray of all the data in
    madrigalFile from a given time

    Inputs:
    madrigalFile - path to Madrigal LOS Hdf5 file
    timeArr - a numpy array of the times (ut1_unix) in the entire
    file.
    unixTime - the unix time to select

    Returns a numpy recarray of the subset of data from the file from
    that one site
    """
    with h5py.File(madrigalFile, "r") as f:
        # get a list of all the indices with the right site
        time_indices = timeArr == unixTime
        timeData = f["Data"]["Table Layout"][time_indices]
        return timeData


p = argparse.ArgumentParser(description="Read Madrigal LOS file")
p.add_argument("filename", help="Madrigal LOS file")
args = p.parse_args()

fn = args["filename"]

t = time.time()

# we only need to do this step once, so that getting data from other sites is faster
with h5py.File(fn, "r") as f:
    timeArr = f["Data"]["Table Layout"]
    # get the list of unique times from this array
    times = numpy.unique(timeArr)

print(f"times in file {fn}:", times)

print(timeArr)

timeData = getTimeData(fn, timeArr, times[0])
print(
    f"Took {time.time() - t} secs to get {len(timeData)} measurements from the first time"
)
# show that accessing the second time is faster
t = time.time()
timeData = getTimeData(fn, timeArr, times[1])
print(
    f"{time.time() - t} secs to get {len(timeData)} measurements from the second time"
)
print("The following are the column names and data types in this file:")
for colName, colType in timeData.dtype.fields.items():
    print(f"{colName}\t{colType}")
print("For example, here are the first four rows:")
for i in range(4):
    print(i, timeData[i])
print("For example, here are the first four rows of line of site TEC:")
for i in range(4):
    print(i, timeData["los_tec"][i])
print(f"There are {len(times)} unique times in the file")