scivision/transcarread

View on GitHub
transcarread/io.py

Summary

Maintainability
A
0 mins
Test Coverage
from pathlib import Path
from typing import Dict, Any, Tuple
from datetime import datetime, timedelta
import numpy as np


def parseionoheader(h: np.ndarray) -> Dict[str, Any]:
    """
    for reading 90kmmaxpt123.dat, which may have been generated by MSIS90
    This file was figured out by inspection of dir.source/transconvec_13.op.f

    Parameters
    ----------

    dz: altitude step

    variables:
    nx: number of altitudes
    raw: 2-D array of data for each altitude, this is what we'll interpolate and
    write back to disk

    Note: by inspection, data values of 90kmmaxpt123.dat start at byte 504

    TRANSCAR takes the initial conditions file as an initial state at t=0, and then
    computes the state of the ionosphere for times t=0+{T1,T2,T3....Tn}.
    This output is stored in dir.output/transcar_output, which is read by
    read_tra.py in a similar manner to this code
    """
    assert 1 <= h[3] <= 12
    assert 1 <= h[4] <= 31
    assert 0 <= h[5] < 24
    assert 0 <= h[6] < 60
    assert 0 <= h[7] < 60
    # ... and so on with asserts. Just checking we aren't reading the totally wrong type of file
    # not a Series because all have to be same datatype
    # not a Dataframe because it's only 1-D
    hd = {
        "nx": h[0].astype(int),
        "ncol": h[1].astype(int),
        "intpas": h[8],
        "longeo": h[9],
        "latgeo": h[10],
        "lonmag": h[11],
        "latmag": h[12],
        "tmag": h[13],
        "f1072": h[14],
        "f1073": h[15],
        "ap2": h[16],
        "ikp": h[17],
        "dTinf": h[18],
        "dUinf": h[19],
        "cofo": h[20],
        "cofh": h[21],
        "cofn": h[22],
        "chi": h[23],
        "approx": h[36],
    }

    # h[37] last non-zero value till h[59], then zeros till start of data at byte 504
    # h[59] has value of 1.0

    hd["htime"] = datetime(year=h[2], month=h[3], day=h[4], hour=h[5], minute=h[6], second=h[7])

    return hd


def readionoheader(tcofn: Path, nhead: int) -> Tuple[Dict[str, Any], np.ndarray]:
    """ reads BINARY transcar_output file """
    tcofn = Path(tcofn).expanduser()  # not dupe, for those importing externally

    if tcofn.is_dir():  # need this on windows
        raise IsADirectoryError(tcofn)

    with tcofn.open("rb") as f:
        h = np.fromfile(f, np.float32, nhead)

    return parseionoheader(h), h


def readTranscarInput(infn: Path) -> Dict[str, Any]:
    """
    The transcar input file is indexed by line number --this is what the Fortran
      #  code of transcar does, and it's what we do here as well.
    """
    infn = Path(infn).expanduser()
    hd: Dict[str, Any] = {}

    with infn.open("r") as f:
        hd["kiappel"] = int(f.readline().split()[0])
        hd["precfile"] = f.readline().split()[0]
        hd["dtsim"] = float(f.readline().split()[0])  # "dto"
        hd["dtfluid"] = float(f.readline().split()[0])  # "sortie"
        hd["iyd_ini"] = int(f.readline().split()[0])
        hd["dayofsim"] = datetime.strptime(str(hd["iyd_ini"]), "%Y%j")
        hd["simstartUTCsec"] = float(f.readline().split()[0])  # "tempsini"
        hd["simlengthsec"] = float(f.readline().split()[0])  # "tempslim"
        hd["jpreci"] = int(f.readline().split()[0])
        # transconvec calls the next two latgeo_ini, longeo_ini
        hd["latgeo_ini"], hd["longeo_ini"] = [float(a) for a in f.readline().split(None)[0].split(",")]
        # from transconvec, time before precip
        hd["tempsconv_1"] = float(f.readline().split()[0])
        # from transconvec, time after precip
        hd["tempsconv"] = float(f.readline().split()[0])
        hd["step"] = float(f.readline().split()[0])
        # transconvec calls this "postinto"
        hd["dtkinetic"] = float(f.readline().split()[0])
        hd["vparaB"] = float(f.readline().split()[0])
        hd["f107ind"] = float(f.readline().split()[0])
        hd["f107avg"] = float(f.readline().split()[0])
        hd["apind"] = float(f.readline().split()[0])
        hd["convecEfieldmVm"] = float(f.readline().split()[0])
        hd["cofo"] = float(f.readline().split()[0])
        hd["cofn2"] = float(f.readline().split()[0])
        hd["cofo2"] = float(f.readline().split()[0])
        hd["cofn"] = float(f.readline().split()[0])
        hd["cofh"] = float(f.readline().split()[0])
        hd["etopflux"] = float(f.readline().split()[0])
        hd["precinfn"] = f.readline().split()[0]
        hd["precint"] = int(f.readline().split()[0])
        hd["precext"] = int(f.readline().split()[0])
        hd["precipstartsec"] = float(f.readline().split()[0])
        hd["precipendsec"] = float(f.readline().split()[0])

        # derived parameters not in datcar file
        hd["tstartSim"] = hd["dayofsim"] + timedelta(seconds=hd["simstartUTCsec"])
        # TODO verify this isn't added to start
        hd["tendSim"] = hd["dayofsim"] + timedelta(seconds=hd["simlengthsec"])
        hd["tstartPrecip"] = hd["dayofsim"] + timedelta(seconds=hd["precipstartsec"])
        hd["tendPrecip"] = hd["dayofsim"] + timedelta(seconds=hd["precipendsec"])

    return hd